 /*
  * Khoros: $Id: ldggauss.c,v 1.3 1992/03/20 23:26:31 dkhoros Exp $
  */

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

 /*
  * $Log: ldggauss.c,v $
 * Revision 1.3  1992/03/20  23:26:31  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 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: ldggauss.c
 >>>>
 >>>>      Program Name: dggauss
 >>>>
 >>>> Date Last Updated: Sat Feb 15 16:29:36 1992 
 >>>>
 >>>>          Routines: ldggauss - the library call for dggauss
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#define KERNSIZE 2
/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldggauss - library call for dggauss
*
* Purpose:
*    
*    Creates a gaussian noise data set.
*    
*    

* Input:
*    
*    image          as an optional input, this argument can be used to
*                   pass  the  pointer  to  the parameter "kernel", if
*                   this argument is NULL, then the  parameter  option
*                   is not used.
*    
*    numvec         number of data vectors
*    
*    numpts         number of points per data vector
*    
*    mean           mean of data set
*    
*    variance       variance of data set
*    
*    type           data type of output image.  This  must  be  either
*                   VFF_TYP_FLOAT or VFF_TYP_COMPLEX.
*    
*    procdir        process direction.  This can be either 0 for  vec-
*                   tor processing or 1 for band processing.
*    
*    

* Output:
*    
*    image          holds the result of the data generation. The  out-
*                   put data type is the same as the input data type.
*    
*    return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Jeremy Worley
*    
*    

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


/* -library_def */
int ldggauss(image,numvec,numpts,mean,variance,type,procdir)
    struct xvimage **image;
    int    numvec,numpts,procdir,type;
    float  mean,variance;
/* -library_def_end */

/* -library_code */
{
    char *program="ldggauss";
  char  *comment="noise data courtesy of ldggauss.";
  int   i,rows,cols,bands,kernargs,kernsize;
  float *kern, **ensemble;
  int   dunload_vector(),gsin();
  struct xvimage *createimage();
/*
** if image is not NULL, then it should contain a kernel that describes
** the ensemble we want to create.  Copy the image data into kern for
** use later.  Then free up the image pointer so it can be recreated
** below.
** get information from kernel if necessary:
*/

  if(*image!=NULL){
     kern = (float *)(*image)->imagedata;
     if(kern==NULL){
        fprintf(stderr,"%s:  Kernel pointer is NULL.\n",program);
        return(0);
     }
     kernsize = (*image)->num_data_bands;
     kernargs = (*image)->row_size * (*image)->col_size;
     if(kernsize<KERNSIZE || kernargs<numvec){
        fprintf(stderr,"ldgsin:  size of kernel is too small.\n");
        return(0);
     }
     freeimage(*image);
  }else{
    kern = NULL;
  }

/*
** Immediately check bounds on arguments. 
*/

  if(numvec<1 || numpts<1){
     fprintf(stderr,"%s:  Image dimension out of bounds.\n",program);
     return(0);
  }else if(variance<0.0 && kern==NULL){
     fprintf(stderr,"%s:  Variance out of bounds.\n",program);
     return(0);
  }else if(type!=VFF_TYP_FLOAT && type!=VFF_TYP_COMPLEX){
     fprintf(stderr,"%s:  Illegal image data type specified.\n",program);
     return(0);
  }

/*
** setup dimensions based on direction of processing
*/

  cols = 1;
  if(procdir==DSP_VECTOR){
     rows = numvec;
     bands = numpts;
  }else{
     rows = numpts;
     bands = numvec;
  }

/*
** Generate image
*/
  
  *image = createimage(rows,cols,type,1,bands,comment,
                      (unsigned long)0, /* map_row_size */
                      (unsigned long)0, /* map_col_size */
                      (unsigned long)VFF_MS_NONE,
                      (unsigned long)VFF_MAPTYP_NONE,
                      (unsigned long)VFF_LOC_IMPLICIT,
                      (unsigned long)0);/* location dimension */

  if(*image==NULL){
     fprintf(stderr,"%s:  Unable to allocate memory for image.\n",program);
     return(0);
  }

/*
** allocate "final resting place" of data
*/

  ensemble = (float **)malloc(numvec*sizeof(float));
  if(ensemble==NULL){
      fprintf(stderr,"%s:  Unable to allocate work memory\n",program);
      freeimage(*image);
      return(0);
  }

  for(i=0;i<numvec;i++){
      ensemble[i] = (float *)malloc(numpts*sizeof(float));
      if(ensemble[i]==NULL){
         fprintf(stderr,"%s:  Unable to allocate work memory\n",program);
         freeimage(*image);
         for(i=0;i<numvec;i++)if(ensemble[i]!=NULL)free(ensemble[i]);
         free(ensemble);
         return(0);
      }
  }

/*
** make call to mid level routine that generates signal ensemble regardles
** of whether or not the kern has information in it.
*/

  if(!gengauss(ensemble,kern,kernargs,numvec,numpts,mean,variance,type)){
     fprintf(stderr,"%s:  Error during generation of signal.\n",program);
     freeimage(*image);
     for(i=0;i<numvec;i++)if(ensemble[i]!=NULL)free(ensemble[i]);
     free(ensemble);
     return(0);
  }

/*
** reorganize data and place it in image
*/
        
  if(!dunload_vector((char **)ensemble,*image,(unsigned long)type,numvec,
        numpts,procdir)){
      fprintf(stderr,"%s:  Unable to unload vector.\n",program);
      freeimage(*image);
      for(i=0;i<numvec;i++)if(ensemble[i]!=NULL)free(ensemble[i]);
      free(ensemble);
      return(0);
  }

/*
** clean up mess
*/

  for(i=0;i<numvec;i++)if(ensemble[i]!=NULL)free(ensemble[i]);
  free(ensemble);

/*
** normal return
*/

  return(1);
} /* end of ldggauss */

/**************************************************************
*
* MODULE NAME: gengauss
*
*     PURPOSE: generates a 1d signal ensemble
*
*       INPUT:
*
*      OUTPUT: does the work
*
* CALLED FROM: ldggauss
*
* ROUTINES CALLED:
*
**************************************************************/

int gengauss(ensemble,kern,kernargs,numvec,numpts,mean,variance,type)
int            numvec,numpts,type,kernargs;
float          **ensemble,*kern,mean,variance;
{
  char *program = "ldggauss.gengauss";
  float *rpart,*ipart;
  int j,i,siz;
  
/*
** malloc appropriate amount of memory for real part of signal.
** if there are any errors along the way, free up whatever memory was 
** allocated before exiting.
*/

  siz=1;
  if(type==VFF_TYP_COMPLEX)siz=2;

  rpart = (float *)malloc(numpts*sizeof(float));
  if(rpart==NULL){
     fprintf(stderr,"%s:  Unable to malloc temporary workspace.\n",program);
     return(0);
  }

   
/*
** if this is a complex signal ensemble, malloc appropriate amount of 
** memory for the imaginary part.  if there are any errors along the way, 
** free up whatever memory was allocated before exiting.
*/
 
  if(type==VFF_TYP_COMPLEX){
     ipart = (float *)malloc(numpts*sizeof(float));
     if(ipart==NULL){
        fprintf(stderr,"%s:  Unable to malloc temporary workspace.\n",program);
        free(rpart);
        return(0);
     }
  }

/*
** get information from kernel as necessary and loop through all elements
** of the kernel.
*/

  if(kern!=NULL){
     for(i=0;i<numvec;i++){
         mean = kern[i];
         variance = kern[kernargs + i];
         
         if(variance<0.0){
            fprintf(stderr,"%s:  Kernel variance out of bounds\n",
                           program);
            return(0);
         }
          
    
/*
** call low level sine wave generator for real part.
*/

         j=ggauss(rpart,numpts,1,mean,variance);
         if(j==0){
            fprintf(stderr,"%s:  Error during signal generation\n",program);
            free(rpart);
            if(type==VFF_TYP_COMPLEX)free(ipart);
            return(0);
         }

/*
** call low level sine wave generator for imaginary part.
*/

         if(type==VFF_TYP_COMPLEX){
            j=ggauss(ipart,numpts,1,mean,variance);
            if(j==0){
               fprintf(stderr,"%s:  Error during signal generation\n",
                              program);
               free(rpart);
               free(ipart);
               return(0);
            }
         }

/*
** interlace real and imaginary parts if necessary
** reorganize data and place it in ensemble
*/

         if(type==VFF_TYP_COMPLEX){
            for(j=0;j<numpts;j++){
                ensemble[i][2*j] = rpart[j];
                ensemble[i][2*j+1] = ipart[j];
            }
         }else{
            bcopy(rpart,ensemble[i],numpts*sizeof(float));
         }
     } /* end for loop */

/*************************************************************************
** Now for the case where kern == NULL
**************************************************************************/

  }else{

/*
** call low level sine wave generator for real part.
*/

     j=ggauss(rpart,numpts,1,mean,variance);
     if(j==0){
        fprintf(stderr,"%s:  Error during signal generation\n",program);
        free(rpart);
        if(type==VFF_TYP_COMPLEX)free(ipart);
        return(0);
     }

/*
** call low level sine wave generator for imaginary part.
*/

     if(type==VFF_TYP_COMPLEX){
        j=ggauss(ipart,numpts,1,mean,variance);
        if(j==0){
           fprintf(stderr,"%s:  Error during signal generation\n",program);
           free(rpart);
           free(ipart);
           return(0);
        }
     }



/*
** interlace real and imaginary parts if necessary
** reorganize data and place it in ensemble
*/

     if(type==VFF_TYP_COMPLEX){
        for(j=0;j<numpts;j++){
            ensemble[0][2*j] = rpart[j];
            ensemble[0][2*j+1] = ipart[j];
        }
     }else{
        bcopy(rpart,ensemble[0],numpts*sizeof(float));
     }

     for(i=1;i<numvec;i++)bcopy((char *)ensemble[0],ensemble[i],
                                siz*numpts*sizeof(float)); 

  } /* end of big if statement */

/*
** clean up mess
*/ 

  free(rpart);
  if(type==VFF_TYP_COMPLEX)free(ipart);

/*
** normal return
*/

  return(1);

}
/* -library_code_end */
