// gauss.c
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.
// This program was copied from "Numerical recipes in C",
//    used to generate Gauss noise signal  */
// This code was originally developed by: Juan Jose Murillo
//    modifications by: Ralph J. Jackson  <rjj@cis.ufl.edu>

#include <stream.h>
#include <math.h>

#include "gauss.h"


#define  M1   259200        /* constant used by ran1() */
#define  IA1  7141
#define  IC1  54773
#define  RM1  (1.0/M1)
#define  M2   134456
#define  IA2  8121
#define  IC2  28411
#define  RM2  (1.0/M2)
#define  M3   243000
#define  IA3  4561
#define  IC3  51349


//float    ran1();         /* uniform noise generation function */

//_______________________________________________________________
// Return a uniform deviate between 0.0 and 1.0. set idum to any
//     negtive value to initialize or reinitialize the sequence
// Copied from "Numerical Recipes in C"

float ran1( int *idum )
{
    static long ix1,ix2,ix3;
    static float r[98];
    float temp;
    static int iff=0;
    int j;
    
    if( *idum<0 || iff == 0 ) {
	iff = 1;
	ix1 = (IC1-(*idum))%M1;
	ix1 = (IA1*ix1+IC1)%M1;
	ix2 = ix1%M2;
	ix1 = (IA1*ix1+IC1)%M1;
	ix3 = ix1%M3;
	for( j=1; j<=97; j++ ) {
	    ix1 = (IA1*ix1+IC1)%M1;
	    ix2 = (IA2*ix2+IC2)%M2;
	    r[j] = (ix1+ix2*RM2)*RM1;
	}
	*idum = 1;
    }
    ix1 = (IA1*ix1+IC1)%M1;
    ix2 = (IA2*ix2+IC2)%M2;
    ix3 = (IA3*ix3+IC3)%M3;
    j = (int)(1+((97*ix3)/M3));
    if(j>97 || j<1)  {
	fprintf(stderr,"RAN1: This can not happen.j= %d\n",j);
	/* exit(1);*/
    }
    temp = r[j];
    r[j] = (ix1+ix2*RM2)*RM1;
    return temp;
}


//_______________________________________________________________
// Returns a normally distributed deviate with zero mean and
//      unit variance, using ran1(idum) as source of uniform deviates
// Copied from "Numerical Recipes in C"

float gasdev( int *idum )
{
    static int iset = 0;
    static float gset;
    float fac,r,v1,v2;

    if ( iset == 0 )  {
	do {
	    v1 = 2.0*ran1(idum)-1.0;
	    v2 = 2.0*ran1(idum)-1.0;
	    r  = v1*v1+v2*v2;
	}
	while(r >= 1.0);

	fac = sqrt(-2.0*log(r)/r);
	gset = v1*fac;
	iset = 1;
	return v2*fac;
    }
    else {
	iset = 0;
	return gset;
    }
}


// The five parameters for the constructor are :
// mean, standard deviation, percentage of the image perturbed by noise,
// min and max values that the image will take.


IA_Gauss::IA_Gauss( float _mean, float _stdv,
		    float _p_noise, float _min, float _max)
:mean(_mean), stdv(_stdv),
 p_noise(_p_noise), min(_min), max(_max) {}


IA_Closure<float, float> *
IA_Gauss::clone_self() const {
    return new IA_Gauss(*this);
}

float
IA_Gauss::operator() (const float &v) const {
    int n;
    if (ran1(&n) < p_noise) {
	float newv = float(stdv * gasdev(&n) + mean + v);
	if (newv > max) newv = max;
	if (newv < min) newv = min;
	return newv;
    }
    else
	return v;
}
