 /*
  * Khoros: $Id: lvgauss.c,v 1.2 1991/12/18 09:23:23 dkhoros Exp $
  */

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

 /*
  * $Log: lvgauss.c,v $
 * Revision 1.2  1991/12/18  09:23:23  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: lvgauss.c
 >>>>
 >>>>      Program Name: vgauss
 >>>>
 >>>> Date Last Updated: Mon Dec  9 23:22:24 1991 
 >>>>
 >>>>          Routines: lvgauss - the library call for vgauss
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


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


/****************************************************************
*
* Routine Name: lvgauss - library call for vgauss
*
* Purpose:
*    
*    Create an image containing one or more Gaussians
*    
*    

* Input:
* Output:
*
* Written By: Scott Wilson
****************************************************************/


/* -library_def */
int
lvgauss(image,r,c,x,y,vx,vy,ra,amp,n,params)
int r,c,n;
float x,y,vx,vy,ra,amp;
struct xvimage *params,**image;
/* -library_def_end */

/* -library_code */
{
  float xloc,yloc;
  float varx,vary;
  float rotang;
  float ampl;
  int i,j,ng,k;
  struct xvimage *img,*createimage();
  float *data,*par;
  float geng();

  /* Create the output image */
  img   = createimage((unsigned long)r,  /* number of rows */
                     (unsigned long)c,   /* number of columns */
                     (unsigned long)VFF_TYP_FLOAT,   /* data storage type */
                     (unsigned long)1,      /* num_of_images */
                     (unsigned long)1,      /* num_data_bands */
                     "created by vgauss",    /* 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 (img == NULL)
    {
      fprintf(stderr,"lvgauss: Unable to allocate new image!\n");
      *image = NULL;
      return(0);
    }


  /* Fill in the Gaussian values */
  data = (float *)(img->imagedata);
  if (params == NULL)
    {
      /* If only want the single Gaussian */
      for (i=0; i<r; i++)
        for (j=0; j<c; j++)
          data[i*c+j] = geng(x,y,vx,vy,ra,amp,n,j,i);
    }
  else
    {
      /* Must use the parameters from the parameter image */
      if (params->num_data_bands != 6)
        {
          fprintf(stderr,"lvgauss: Parameter image must have 6 bands\n");
          freeimage(img);
          return(0);
        }
      ng = params->row_size*params->col_size;
      par = (float *)(params->imagedata);
      bzero(data,r*c*sizeof(float));
      for (k=0; k<ng; k++)
        {
          xloc   = par[0*ng+k];
          yloc   = par[1*ng+k];
          varx   = par[2*ng+k];
          vary   = par[3*ng+k];
          rotang = par[4*ng+k];
          ampl   = par[5*ng+k];
          for (i=0; i<r; i++)
            for (j=0; j<c; j++)
              data[i*c+j] += geng(xloc,yloc,varx,vary,rotang,ampl,n,j,i);
        }
    }
  *image = img;
  return(1);
}

#define DEG2RAD 1.74533e-2
#define PI 3.1415926

float
geng(xloc,yloc,varx,vary,rotang,ampl,n,x,y)
float xloc,yloc,varx,vary,rotang,ampl;
int n,x,y;
 {
   /* Compute the value of the Gaussian distribution at the (x,y) location
      given using the supplied parameters for the distribution */
   float p,sp;     /* Correlation and sqrt(1-p^2) */
   float s1,s2;    /* Standard deviations */
   float a1,a2,a3; /* x^2, cross, and y^2 terms */
   float a4;       /* Exponent multiplier */
   float a5;       /* Full exponent */
   float a6;       /* Constant in front */
   float f;        /* Result */

   if (n)
     {
       /* Want normalized Gaussian */
       p = sin((double)(DEG2RAD*rotang));
       sp = sqrt((double)(1.0-p*p));
       s1 = sqrt((double)varx);
       s2 = sqrt((double)vary);

       a1 = (x-xloc)*(x-xloc)/varx;
       a2 = 2.0*p*(x-xloc)*(y-yloc)/(s1*s2);
       a3 = (y-yloc)*(y-yloc)/vary;
       a4 = 1.0/(2.0*(1.0-p*p));
       a5 = -a4*(a1-a2+a3);
       a6 = 1.0/(2.0*PI*s1*s2*sp);
       f = a6*exp((double)a5);
       return(f);
     }
   else
     {
       /* Use the given amplitude  (un-normalized) */
       p = sin((double)(DEG2RAD*rotang));
       s1 = sqrt((double)varx);
       s2 = sqrt((double)vary);

       a1 = (x-xloc)*(x-xloc)/varx;
       a2 = 2.0*p*(x-xloc)*(y-yloc)/(s1*s2);
       a3 = (y-yloc)*(y-yloc)/vary;
       a4 = 1.0/(2.0*(1.0-p*p));
       a5 = -a4*(a1-a2+a3);
       f = ampl*exp((double)a5);
       return(f);
     }
 }
/* -library_code_end */
