 /*
  * Khoros: $Id: lvgfractal.c,v 1.2 1991/12/18 09:24:02 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvgfractal.c,v 1.2 1991/12/18 09:24:02 dkhoros Exp $";
#endif

 /*
  * $Log: lvgfractal.c,v $
 * Revision 1.2  1991/12/18  09:24:02  dkhoros
 * HellPatch3
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * 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 to 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: lvgfractal.c
 >>>>
 >>>>      Program Name: vgfractal
 >>>>
 >>>> Date Last Updated: Mon Dec  9 23:26:43 1991 
 >>>>
 >>>>          Routines: lvgfractal - the library call for vgfractal
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#define TWOPI 6.283185307
double uno;
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvgfractal - library call for vgfractal
*
* Purpose:
*    
*    Creates a fractal  image  using  the  midpoint  displacement
*    algorithm.
*    
*    

* Input:
*    
*    image          xvimage structure
*    
*    maxlevel       maximum number of recursions
*    
*    hurst          desired hurst coeffient of image
*    
*    sigma          standard deviation of random numbers
*    
*    seed           seed for random number generator
*    
*    add_flag       flag to turn random additions on/off
*    
*    

* Output:
*    
*    image          output xvimage structure. Holds the resulting  the
*                   fractal image.
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Charlie Gage
*    
*    

****************************************************************/


/* -library_def */
int
lvgfractal(image, maxlevel, hurst, sigma, seed, coast, add_flag)
struct xvimage **image;
int maxlevel, seed, coast, add_flag;
float  hurst, sigma;
/* -library_def_end */

/* -library_code */
{
    int     i,                      /* general indexing variables */
            j,
            d,
            x,                      /* index for cols */
            y,                      /* index for rows */
            rows,                   /* number of rows in image */
            cols,                   /* number of cols in image */
            s,                      /* step size variable */
            n,                      /* array size, from # of recursions */
            stage;                  /* stage index for recursions */ 

    float   delta,                  /* standard deviation of random values */ 
            H,                      /* hurst coefficient */
            x0,                     /* temporary variables */
            x1,
            x2,
            x3,
            gauss(),                /* function definition */
            *fptr,                  /* pointer to data area*/
            **result;               /* 2-D array for new image data */

    double  num_recur,              /* number of recursions */
            max_recur;

    struct xvimage *kern, *createimage();

    uno = 2.0 * exp(0.825);
                            /* initialize random number generator */
    if (seed != 0)
    {
      (void) vsrandom(seed);
    }
    else
    {
      (void) vsrandom((int)time(NULL));
    }

    delta = sigma;
    H = hurst;
    max_recur = (double) maxlevel;
    num_recur = pow(2.0, max_recur);  /* get max number of recursions */
    n = (int) num_recur;
    rows = n + 1;                     /* determine row size */
    cols = n + 1;                     /* determine col size */


                  /* Call createimage to set up space for a new image */
    kern = createimage((unsigned long) rows,    /* number of rows */
           (unsigned long) cols,   /* number of columns */
           (unsigned long)         /* data_storage_type */
           VFF_TYP_FLOAT,
           (unsigned long) 1,   /* num_of_images */ 
                        (unsigned long) 1,      /* num_data_bands */ 
           "created by vgfractal",      /* comment */
                        (unsigned long) 0,      /* map_row_size */
                        (unsigned long) 0,      /* map_col_size */
                        (unsigned long)
                        VFF_MS_NONE,             /* map_scheme */
                        (unsigned long)
                        VFF_MAPTYP_NONE,        /* map_storage_type */
                        (unsigned long)
                        VFF_LOC_IMPLICIT,       /* location_type */
                        (unsigned long) 0);     /* location_dim */

    if (kern == NULL)
    {
       (void)fprintf(stderr,"lvgfractal: Unable to allocate new image!\n");
       *image = NULL;
       return(0);
    }

    result = (float **) malloc((unsigned int)cols * rows * sizeof(float));
    if (result == NULL)
    {
     (void)fprintf(stderr,"lvgfractal: Unable to allocate space for new image!\n");
     return(0);
    }

    for (i = 0; i < rows; i++)
    {
      result[i] = (float *) malloc((unsigned int)cols * sizeof(float));
      if (result[i] == NULL)
      {
       (void)fprintf(stderr,"lvgfractal: Unable to allocate space for new image!\n");
       return(0);
      }
    }

    fptr = (float *) kern->imagedata;       /* set pointer to data area */

   /*-----------------------------------------------------------------*
    * generate the fractal image using midpoint displacement          *
    *-----------------------------------------------------------------*/

                           /* set corners to initial random values */
    result[0][0] = delta * gauss();
    result[0][n] = delta * gauss();
    result[n][0] = delta * gauss();
    result[n][n] = delta * gauss();
    s = n;
    d = n / 2;

    for (stage = 0; stage < maxlevel; stage++)
    {
      delta = delta * pow(0.5, 0.5*H);
      for (x = d; x <= (n-d);) 
      {
        for (y = d; y <= (n-d);) 
        {
          x0 = result[x+d][y+d];
          x1 = result[x+d][y-d];
          x2 = result[x-d][y+d];
          x3 = result[x-d][y-d];
          result[x][y] = (x0 + x1 + x2 + x3)/ 4.0 + delta * gauss();
          y = y + s;         /* step s */
        }
        x = x + s;           /* step s */
      }
                            /* displace other points if needed */
    if (add_flag)
    {
      for (x = 0; x <= n;)
      {
        for (y = 0; y <= n;)
        {
          result[x][y] = result[x][y] + delta * gauss();
          y = y + s;        /* step s */
        }
        x = x + s;          /* step s */
      }
    }
                            /* change from grid type II to type I */
    delta = delta * pow(0.5, 0.5*H);
                            /* interpolate and offset boundary grid points */
    for (x = d; x <= (n-d);)
    {
      x0 = result[x+d][0];
      x1 = result[x-d][0];
      x2 = result[x][d];
      result[x][0] = (x0 + x1 + x2)/3.0 + delta * gauss();
      x0 = result[x+d][n];
      x1 = result[x-d][n];
      x2 = result[x][n-d];
      result[x][n] = (x0 + x1 + x2)/3.0 + delta * gauss();
      x0 = result[0][x+d];
      x1 = result[0][x-d];
      x2 = result[d][x];
      result[0][x] = (x0 + x1 + x2)/3.0 + delta * gauss();
      x0 = result[n][x+d];
      x1 = result[n][x-d];
      x2 = result[n-d][x];
      result[n][x] = (x0 + x1 + x2)/3.0 + delta * gauss();
      x = x + s;     /* step s */
    }
                      /* interpolate and offset interior grid points */
    for (x = d; x <= (n-d);)
    {
       for (y = s; y <= (n-d);) 
       {
         x0 = result[x][y+d];
         x1 = result[x][y-d];
         x2 = result[x+d][y];
         x3 = result[x-d][y];
         result[x][y] = (x0 + x1 + x2 + x3)/ 4.0 + delta * gauss();
         y = y + s;         /* step s */
       }
       x = x + s;           /* step s */
    }
    for (x = s; x <= (n-d);)
    {
       for (y = d; y <= (n-d);) 
       {
         x0 = result[x][y+d];
         x1 = result[x][y-d];
         x2 = result[x+d][y];
         x3 = result[x-d][y];
         result[x][y] = (x0 + x1 + x2 + x3)/ 4.0 + delta * gauss();
         y = y + s;         /* step s */
       }
       x = x + s;           /* step s */
    }
                          /* displace other points if needed */
    if (add_flag)
    {
      for (x = 0; x <= n;)
      {
        for (y = 0; y <= n;)
        {
          result[x][y] = result[x][y] + delta * gauss();
          y = y + s;        /* step s */
        }
        x = x + s;          /* step s */
      }
      for (x = d; x <= (n-d);)
      {
        for (y = d; y <= (n-d);)
        {
          result[x][y] = result[x][y] + delta * gauss();
          y = y + s;        /* step s */
        }
        x = x + s;          /* step s */
      }
    }

   s = s / 2;
   d = d / 2;

 }    /* end of for(stage) */

               /* put data in a 1-D VIFF image array format */
 for (i = 0; i < rows; i++)
 {
   for (j = 0; j < cols; j++)
   {
     if ((coast) && (result[i][j] < 0.0))
       fptr[i * cols + j] = 0.0;
     else 
       fptr[i * cols + j] = result[i][j];
   }
 }

 *image = kern;
 return(1);

}  /* End of lvgfractal */


/*-------------------------------------------------------------*
 *  routine to generate a gaussian random number
 *-------------------------------------------------------------*/

float gauss()
{
  long    big, u1, u2, vrandom();
  double  temp, rn, var, mean;
  float   f, f1, f2;

  var = 1.0;              /* unit variance and zero mean */
  mean = 0.0;

  u1 = vrandom(&big);
  u2 = vrandom(&big);
  f1 = (float) u1/big;
  f2 = (float) u2/big;
  temp = sqrt((double)(-uno*var*log10((double)f1)));
  rn = temp * cos(TWOPI*f2) + mean;

  f = (float)rn;

  return(f);
}

/* -library_code_end */
