Article 13643 of alt.sources:

Archive-name: pnmaddnoise.c_1.0
Submitted-by: burns@cac.psu.edu

pnmaddnoise is a pbmplus/netpbm program that adds one of several types
of noise to an image (impulse, gaussian, etc).  Included in the shar 
file is the source code, man page and README file for pnmaddnoise.

- Mike

----------------------------------------------------------------------------
Mike Burns                                                UNIX Systems Group 
burns@cac.psu.edu                              Center for Academic Computing 
(814) 863-5606                             The Pennsylvania State University

---- Cut Here and unpack ----
#!/bin/sh
# shar:	Shell Archiver  (v1.22)
#
#	Run the following text with /bin/sh to create:
#	  pnmaddnoise.README
#	  pnmaddnoise.1
#	  pnmaddnoise.c
#
sed 's/^X//' << 'SHAR_EOF' > pnmaddnoise.README &&
XREADME for pnmaddnoise 1.0
X
Xpnmaddnoise is a pbmplus/netpbm program that adds one of the following types 
Xof noise to an image.
X
X	gaussian
X	multiplicative gaussian
X	impulse
X	laplacian
X	poisson
X
XThis was written partly to generate noisy images to test the median filter 
XI was working on.  I found many descriptions of noise, but not too many
Xdescriptions of how to add noise to an image, so it is quite possible this
Xdoesn't implement the noise in the best fashion, but the resultant images
Xcertainly do have noise in them.
X
XInstallation
X------------
XCopy the files to your {pbmplus|netpbm}/pnm directory and modify the 
XImakefile to include pnmaddnoise in MAN1, SRCS, OBJS, BINS and add the
Xtarget line
X
XNormalPbmplusMathProgramTarget(pnmaddnoise)
X
XThen run xmkmf and make.
X
X- Mike
X
X----------------------------------------------------------------------------
XMike Burns                                                UNIX Systems Group 
Xburns@cac.psu.edu                              Center for Academic Computing 
X(814) 863-5606                             The Pennsylvania State University
X
X
SHAR_EOF
chmod 0644 pnmaddnoise.README || echo "restore of pnmaddnoise.README fails"
sed 's/^X//' << 'SHAR_EOF' > pnmaddnoise.1 &&
X.TH pnmaddnoise 1 "14 November 1995"
X.IX pnmaddnoise
X.SH NAME
Xpnmaddnoise - add noise to a portable anymap
X.SH SYNOPSIS
X.B pnmaddnoise
X.RB [ -type
X.IR noise_type ]
X.RB [ -seed
X.IR int ]
X.br
X.B pnmaddnoise
X.B -type 
X.B gaussian
X.RB [ -sigma1
X.IR value ]
X.RB [ -sigma2
X.IR value ]
X.if n .ti +1.2i
X.RB [ -seed
X.IR int ]
X.RI [ pnmfile ]
X.br
X.B pnmaddnoise
X.B -type 
X.B multiplicative_gaussian
X.RB [ -msigma
X.IR value ]
X.if n .ti +1.2i
X.RB [ -seed
X.IR int ]
X.RI [ pnmfile ]
X.br
X.B pnmaddnoise
X.B -type 
X.B impulse
X.RB [ -tolerence
X.RI ratio ]
X.RB [ -seed
X.IR int ]
X.if n .ti +1.2i
X.RI [ pnmfile ]
X.br
X.B pnmaddnoise
X.B -type 
X.B laplacian
X.RB [ -lsigma
X.IR value ]
X.RB [ -seed
X.IR int ]
X.if n .ti +1.2i
X.RI [ pnmfile ]
X.br
X.B pnmaddnoise
X.B -type 
X.B poisson
X.RB [ -lambda
X.IR value ]
X.RB [ -seed
X.IR int ]
X.if n .ti +1.2i
X.RI [ pnmfile ]
X.br
X.SH DESCRIPTION
XAdds the specified noise type to a portable anymap.
XPPM files are treated one color plane at a time, not as a colorspace.
XSee source code for meaning of options particular to each noise type.
X.SH OPTIONS
X.TP
X.B \-lambda value
XUsed for poisson noise only.  The default value is
X.IR 0.05.
X.TP
X.B -lsigma value
XUsed for laplacian noise only.  The default value is
X.IR 10.0.
X.TP
X.B -msigma value
XUsed for mutliplicative gaussian noise only.  The default value is
X.IR 0.5.
X.TP
X.B -seed int
XUsed for all noise types.  Set the initial random number generator seed value.
X.TP
X.B -sigma1 value
XUsed for gaussian noise only.  The default value is
X.IR 4.0.
X.TP
X.B -sigma2 value
XUsed for gaussian noise only.  The default value is
X.IR 20.0.
X.TP
X.B -tolerence ratio
XUsed for impulse noise only.  
XThe default value is
X.IR 0.10. 
XThis means that 5% of all pixel values will be set to 0 and 5% will be set
Xto MAXVAL.
X.TP
X.B -type noise_type
XSelect type of noise by name.
XThe following noise types are available:
Xgaussian, multiplicative_gaussian, impulse, laplacian, poisson.
XOnly enough letters to be unique are required for the noise type option.
XThe default noise type is 
X.IR gaussian.
X.SH "BUGS"
XPPM files are treated one color plane at a time, not as a color space.
X.SH "REFERENCES"
X"Adaptive Image Restoration in Signal-Dependent Noise" by R. Kasturi
XInstitute for Electronic Science, Texas Tech University, 1982
X.TP
X"Digital Image Processing Algorithms" by Ioannis Pitas,
XPrentice Hall, 1993  ISBN 0-13-145814-0
X.SH "SEE ALSO"
Xpgmnoise(1), pnm(5), pnmmedian(1), pnmsmooth(1)
X.SH AUTHOR
XCopyright (C) 1995 by Mike Burns <burns@cac.psu.edu>
X.\" Permission to use, copy, modify, and distribute this software and its
X.\" documentation for any purpose and without fee is hereby granted, provided
X.\" that the above copyright notice appear in all copies and that both that
X.\" copyright notice and this permission notice appear in supporting
X.\" documentation.  This software is provided "as is" without express or
X.\" implied warranty.
SHAR_EOF
chmod 0644 pnmaddnoise.1 || echo "restore of pnmaddnoise.1 fails"
sed 's/^X//' << 'SHAR_EOF' > pnmaddnoise.c &&
X/* pnmaddnoise.c
X** 
X** Add gaussian, multiplicative gaussian, impulse, laplacian or 
X** poisson noise to a portable anymap.
X** 
X** Version 1.0  November 1995
X**
X** Copyright (C) 1995 by Mike Burns (burns@cac.psu.edu)
X**
X** Permission to use, copy, modify, and distribute this software and its
X** documentation for any purpose and without fee is hereby granted, provided
X** that the above copyright notice appear in all copies and that both that
X** copyright notice and this permission notice appear in supporting
X** documentation.  This software is provided "as is" without express or
X** implied warranty.
X*/
X
X/* References
X** ----------
X** "Adaptive Image Restoration in Signal-Dependent Noise" by R. Kasturi
X** Institute for Electronic Science, Texas Tech University  1982
X**
X** "Digital Image Processing Algorithms" by Ioannis Pitas
X** Prentice Hall, 1993  ISBN 0-13-145814-0
X*/
X
X#include "pnm.h"
X#include <math.h>
X
X#ifndef M_PI
X#define M_PI    3.14159265358979323846
X#endif
X
X#define RANDOM_MASK 0x7FFF	/* only compare lower 15 bits.  Stupid PCs. */
X
X#define GAUSSIAN		0
X#define IMPULSE 		1	/* aka salt and pepper noise */
X#define LAPLACIAN		2
X#define MULTIPLICATIVE_GAUSSIAN	3
X#define POISSON			4
X#define MAX_NOISE_TYPES		5
X
X#define CHECK_GRAY \
X    if ( gtemp < 0.0 ) g = 0; \
X    else if ( gtemp > maxval ) g = maxval; \
X    else g = gtemp + 0.5;
X
X#define CHECK_RED \
X    if ( rtemp < 0.0 ) r = 0; \
X    else if ( rtemp > maxval ) r = maxval; \
X    else r = rtemp + 0.5;
X
X#define CHECK_GREEN \
X    if ( gtemp < 0.0 ) g = 0; \
X    else if ( gtemp > maxval ) g = maxval; \
X    else g = gtemp + 0.5;
X
X#define CHECK_BLUE \
X    if ( btemp < 0.0 ) b = 0; \
X    else if ( btemp > maxval ) b = maxval; \
X    else b = btemp + 0.5;
X
Xstatic void gaussian_noise ARGS(( FILE *ifp, float sigma1, float sigma2 ));
Xstatic void impulse_noise ARGS(( FILE *ifp, float tolerence ));
Xstatic void laplacian_noise ARGS(( FILE *ifp, float lsigma ));
Xstatic void multiplicative_gaussian_noise ARGS(( FILE *ifp, float mgsigma ));
Xstatic void poisson_noise ARGS (( FILE *ifp, float lambda ));
X
Xstatic double arand = 32767.0;		/* 2^15-1 in case stoopid computer */
X
Xint 
Xmain( argc, argv )
X    int argc;
X    char *argv[];
X    {
X    FILE *ifp;
X    int argn;
X    int noise_type;
X    int seed;
X    int i;
X    char *usage = "[-type noise_type] [-lsigma x] [-mgsigma x] [-sigma1 x] [-sigma2 x] [-lambda x] [-seed n] [-tolerence ratio] [pgmfile]";
X
X    char *noise_name[] =
X	{ 
X	"gaussian",
X	"impulse",
X	"laplacian",
X	"multiplicative_gaussian",
X	"poisson"
X	};
X    int noise_id[] = 
X	{ 
X	GAUSSIAN,
X	IMPULSE,
X	LAPLACIAN,
X	MULTIPLICATIVE_GAUSSIAN,
X	POISSON
X	};
X    /* minimum number of characters to match noise name for pm_keymatch() */
X    int noise_compare[] =
X	{
X	1,
X	1,
X	1,
X	1,
X	1
X	};
X
X    /* define default values for configurable options */
X    float lambda = 0.05;		
X    float lsigma = 10.0;
X    float mgsigma = 0.5;
X    float sigma1 = 4.0;
X    float sigma2 = 20.0;
X    float tolerence = 0.10;
X
X    pnm_init( &argc, argv );
X
X    seed = time( 0 ) ^ getpid( );
X    noise_type = GAUSSIAN;
X
X    argn = 1;
X    while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
X	{
X	if ( pm_keymatch( argv[argn], "-lambda", 3 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( 
X		    "incorrect number of arguments for -lambda option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -lambda option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    lambda = atof( argv[argn] );
X	    }
X	else if ( pm_keymatch( argv[argn], "-lsigma", 3 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( 
X		    "incorrect number of arguments for -lsigma option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -lsigma option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    lsigma = atof( argv[argn] );
X	    }
X	else if ( pm_keymatch( argv[argn], "-mgsigma", 2 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( 
X		    "incorrect number of arguments for -mgsigma option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -mgsigma option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    mgsigma = atof( argv[argn] );
X	    }
X	else if ( pm_keymatch( argv[argn], "-seed", 3 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( "incorrect number of arguments for -seed option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -seed option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    seed = atoi(argv[argn]);
X	    }
X	else if ( pm_keymatch( argv[argn], "-sigma1", 7 ) ||
X		pm_keymatch( argv[argn], "-s1", 3 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( 
X		    "incorrect number of arguments for -sigma1 option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -sigma1 option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    sigma1 = atof( argv[argn] );
X	    }
X	else if ( pm_keymatch( argv[argn], "-sigma2", 7 ) ||
X		pm_keymatch( argv[argn], "-s2", 3 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( 
X		    "incorrect number of arguments for -sigma2 option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -sigma2 option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    sigma2 = atof( argv[argn] );
X	    }
X	else if ( pm_keymatch( argv[argn], "-tolerence", 3 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( 
X		    "incorrect number of arguments for -tolerence option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -tolerence option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    tolerence = atof( argv[argn] );
X	    }
X	else if ( pm_keymatch( argv[argn], "-type", 3 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( "incorrect number of arguments for -type option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -type option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    /* search through list of valid noise types and compare */
X	    i = 0;
X	    while ( ( i < MAX_NOISE_TYPES ) && !pm_keymatch( argv[argn], 
X		noise_name[i], noise_compare[i] ) )
X		++i;
X	    if ( i >= MAX_NOISE_TYPES )
X		{
X		pm_message( "invalid argument to -type option: %s", 
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    noise_type = noise_id[i];
X	    }
X	else
X	    pm_usage( usage );
X	++argn;
X	}
X
X    if ( argn < argc )
X	{
X	ifp = pm_openr( argv[argn] );
X	argn++;
X	}
X    else
X	ifp = stdin;
X
X    if ( argn != argc )
X	pm_usage( usage );
X
X    srandom( seed );
X
X    switch ( noise_type )
X	{
X	case GAUSSIAN:
X	gaussian_noise( ifp, sigma1, sigma2 );
X	break;
X	
X	case IMPULSE:
X	impulse_noise( ifp, tolerence );
X	break;
X
X	case LAPLACIAN:
X	laplacian_noise( ifp, lsigma );
X	break;
X
X	case MULTIPLICATIVE_GAUSSIAN:
X	multiplicative_gaussian_noise( ifp, mgsigma );
X	break;
X	
X	case POISSON:
X	poisson_noise( ifp, lambda );
X	break;
X
X	default:  
X	pm_error( 
X	    "This should not be happening to you... no noise of that type");
X	break;
X	}
X
X    exit( 0 );
X}
X
X
X
X/* Gaussian Noise
X**
X** Based on Kasturi/Algorithms of the ACM
X*/
X
Xstatic void
Xgaussian_noise( FILE *ifp, float sigma1, float sigma2 )
X    {
X    register xel *xP;
X    register int col;
X    xel *xelrow;
X    xelval r, g, b;
X    xelval maxval;
X    int rows, cols, format, row;
X    double x1, x2, xn, yn;
X    double rtemp, gtemp, btemp;
X
X    int forceplain = 0;
X
X    pnm_readpnminit( ifp, &cols, &rows, &maxval, &format );
X    xelrow = pnm_allocrow( cols );
X
X    pnm_writepnminit( stdout, cols, rows, maxval, format, forceplain );
X
X    for ( row = 0; row < rows; ++row )
X	{
X	pnm_readpnmrow( ifp, xelrow, cols, maxval, format );
X	for ( col = 0, xP = xelrow; col < cols; ++col, ++xP )
X	    {
X	    switch ( PNM_FORMAT_TYPE(format) )
X		{
X		case PPM_TYPE:
X		x1 = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( x1 == 0.0 )
X		    x1 = 1.0;
X		x2 = ( random( ) & RANDOM_MASK ) / arand;
X		xn = sqrt( -2.0 * log( x1 ) ) * cos( 2.0 * M_PI * x2 );
X		yn = sqrt( -2.0 * log( x1 ) ) * sin( 2.0 * M_PI * x2 );
X		r = PPM_GETR( *xP );
X		rtemp = r + ( sqrt( (double) r ) * sigma1 * xn ) + ( sigma2 * yn );
X
X		x1 = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( x1 == 0.0 )
X		    x1 = 1.0;
X		x2 = ( random( ) & RANDOM_MASK ) / arand;
X		xn = sqrt( -2.0 * log( x1 ) ) * cos( 2.0 * M_PI * x2 );
X		yn = sqrt( -2.0 * log( x1 ) ) * sin( 2.0 * M_PI * x2 );
X		g = PPM_GETG( *xP );
X		gtemp = g + ( sqrt( (double) g ) * sigma1 * xn ) + ( sigma2 * yn );
X
X		x1 = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( x1 == 0.0 )
X		    x1 = 1.0;
X		x2 = ( random( ) & RANDOM_MASK ) / arand;
X		xn = sqrt( -2.0 * log( x1 ) ) * cos( 2.0 * M_PI * x2 );
X		yn = sqrt( -2.0 * log( x1 ) ) * sin( 2.0 * M_PI * x2 );
X		b = PPM_GETB( *xP );
X		btemp = b + ( sqrt( (double) b ) * sigma1 * xn ) + ( sigma2 * yn );
X
X		CHECK_GREEN;
X		CHECK_BLUE;
X		PPM_ASSIGN( *xP, r, g, b );
X		break;
X
X
X		default: /* PGM */
X		/* Need to have a non-zero x1, so if random() returns 0, set
X		** it to 1.0 (the maximum value x1 can take).
X		*/
X		x1 = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( x1 == 0.0 )
X		    x1 = 1.0;
X		x2 = ( random( ) & RANDOM_MASK ) / arand;
X		xn = sqrt( -2.0 * log( x1 ) ) * cos( 2.0 * M_PI * x2 );
X		yn = sqrt( -2.0 * log( x1 ) ) * sin( 2.0 * M_PI * x2 );
X		g = PNM_GET1( *xP );
X		gtemp = g + ( sqrt( (double) g ) * sigma1 * xn ) + ( sigma2 * yn );
X
X		CHECK_GRAY;
X		PNM_ASSIGN1( *xP, g );
X		break;
X		}
X	    }
X	pnm_writepnmrow( stdout, xelrow, cols, maxval, format, forceplain );
X	}
X    pnm_freerow( xelrow );
X    }
X
X
X
X/* Impulse (Salt and Pepper) Noise
X**
X*/
X
Xstatic void
Ximpulse_noise( FILE *ifp, float tolerence )
X    {
X    register xel *xP;
X    register int col;
X    xel *xelrow;
X    xelval r, g, b;
X    xelval maxval;
X    int rows, cols, format, row;
X    double rtemp, gtemp, btemp;
X    double sap;
X    float low_tol, high_tol;
X
X    int forceplain = 0;
X
X    pnm_readpnminit( ifp, &cols, &rows, &maxval, &format );
X    xelrow = pnm_allocrow( cols );
X
X    pnm_writepnminit( stdout, cols, rows, maxval, format, forceplain );
X
X    low_tol = tolerence / 2.0;
X    high_tol = 1.0 - ( tolerence / 2.0 );
X
X    for ( row = 0; row < rows; ++row )
X	{
X	pnm_readpnmrow( ifp, xelrow, cols, maxval, format );
X	for ( col = 0, xP = xelrow; col < cols; ++col, ++xP )
X	    {
X	    switch ( PNM_FORMAT_TYPE( format ) )
X		{
X		case PPM_TYPE:
X		sap = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( sap < low_tol ) 
X		    r = 0;
X		else if ( sap >= high_tol )
X		    r = maxval;
X		else
X		    r = PPM_GETR( *xP );
X
X		sap = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( sap < low_tol ) 
X		    g = 0;
X		else if ( sap >= high_tol )
X		    g = maxval;
X		else
X		    g = PPM_GETG( *xP );
X
X		sap = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( sap < low_tol ) 
X		    b = 0;
X		else if ( sap >= high_tol )
X		    b = maxval;
X		else
X		    b = PPM_GETB( *xP );
X
X		PPM_ASSIGN( *xP, r, g, b );
X		break;
X
X		default:
X		sap = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( sap < low_tol ) 
X		    g = 0;
X		else if ( sap >= high_tol )
X		    g = maxval;
X		else
X		    g = PNM_GET1( *xP );
X
X		PNM_ASSIGN1( *xP, g );
X		break;
X		}
X	    }
X	pnm_writepnmrow( stdout, xelrow, cols, maxval, format, forceplain );
X	}
X    pnm_freerow( xelrow );
X    }
X
X
X/* Laplacian Noise 
X**
X** From Pitas' book.
X*/
X
Xstatic void
Xlaplacian_noise( FILE *ifp, float lsigma )
X    {
X    register xel *xP;
X    register int col;
X    xel *xelrow;
X    xelval r, g, b;
X    xelval maxval;
X    int rows, cols, format, row;
X    double rtemp, gtemp, btemp;
X    double u, u1, ime;
X    double INFINITY;
X
X    double EPSILON = 1.0e-5;
X    int forceplain = 0;
X
X    pnm_readpnminit( ifp, &cols, &rows, &maxval, &format );
X    xelrow = pnm_allocrow( cols );
X    INFINITY = (double) maxval;
X
X    pnm_writepnminit( stdout, cols, rows, maxval, format, forceplain );
X
X    for ( row = 0; row < rows; ++row )
X	{
X	pnm_readpnmrow( ifp, xelrow, cols, maxval, format );
X	for ( col = 0, xP = xelrow; col < cols; ++col, ++xP )
X	    {
X	    switch ( PNM_FORMAT_TYPE( format ) )
X		{
X		case PPM_TYPE:
X		u = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( u <= 0.5 )
X		    {
X		    if ( u <= EPSILON )
X			rtemp = PPM_GETR( *xP ) - INFINITY;
X		    else
X			rtemp = PPM_GETR( *xP ) + lsigma * log( 2.0 * u );
X		    }
X		else
X		    {
X		    if ( ( u1 = 1.0 - u ) <= 0.5 * EPSILON )
X			rtemp = PPM_GETR( *xP ) + INFINITY;
X		    else
X			rtemp = PPM_GETR( *xP ) - lsigma * log( 2.0 * u1 );
X		    }
X
X		u = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( u <= 0.5 )
X		    {
X		    if ( u <= EPSILON )
X			gtemp = PPM_GETG( *xP ) - INFINITY;
X		    else
X			gtemp = PPM_GETG( *xP ) + lsigma * log( 2.0 * u );
X		    }
X		else
X		    {
X		    if ( ( u1 = 1.0 - u ) <= 0.5 * EPSILON )
X			gtemp = PPM_GETG( *xP ) + INFINITY;
X		    else
X			gtemp = PPM_GETG( *xP ) - lsigma * log( 2.0 * u1 );
X		    }
X
X		u = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( u <= 0.5 )
X		    {
X		    if ( u <= EPSILON )
X			btemp = PPM_GETB( *xP ) - INFINITY;
X		    else
X			btemp = PPM_GETB( *xP ) + lsigma * log( 2.0 * u );
X		    }
X		else
X		    {
X		    if ( ( u1 = 1.0 - u ) <= 0.5 * EPSILON )
X			btemp = PPM_GETB( *xP ) + INFINITY;
X		    else
X			btemp = PPM_GETB( *xP ) - lsigma * log( 2.0 * u1 );
X		    }
X
X		CHECK_RED;
X		CHECK_GREEN;
X		CHECK_BLUE;
X		PPM_ASSIGN( *xP, r, g, b );
X		break;
X
X		default:
X		u = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( u <= 0.5 )
X		    {
X		    if ( u <= EPSILON )
X			gtemp = PNM_GET1( *xP ) - INFINITY;
X		    else
X			gtemp = PNM_GET1( *xP ) + lsigma * log( 2.0 * u );
X		    }
X		else
X		    {
X		    if ( ( u1 = 1.0 - u ) <= 0.5 * EPSILON )
X			gtemp = PNM_GET1( *xP ) + INFINITY;
X		    else
X			gtemp = PNM_GET1( *xP ) - lsigma * log( 2.0 * u1 );
X		    }
X
X		CHECK_GRAY;
X		PNM_ASSIGN1( *xP, g );
X		break;
X		}
X	    }
X	pnm_writepnmrow( stdout, xelrow, cols, maxval, format, forceplain );
X	}
X    pnm_freerow( xelrow );
X    }
X
X
X/* Multiplicative Gaussian Noise 
X**
X** From Pitas' book. 
X*/
X
Xstatic void
Xmultiplicative_gaussian_noise( FILE *ifp, float mgsigma )
X    {
X    register xel *xP;
X    register int col;
X    xel *xelrow;
X    xelval r, g, b;
X    xelval maxval;
X    int rows, cols, format, row;
X    double rtemp, gtemp, btemp;
X    double uniform, rayleigh, gauss;
X    double INFINITY;
X
X    double EPSILON = 1.0e-5;
X    int forceplain = 0;
X
X    pnm_readpnminit( ifp, &cols, &rows, &maxval, &format );
X    xelrow = pnm_allocrow( cols );
X    INFINITY = (double) maxval;
X
X    pnm_writepnminit( stdout, cols, rows, maxval, format, forceplain );
X
X    for ( row = 0; row < rows; ++row )
X	{
X	pnm_readpnmrow( ifp, xelrow, cols, maxval, format );
X	for ( col = 0, xP = xelrow; col < cols; ++col, ++xP )
X	    {
X	    switch ( PNM_FORMAT_TYPE( format ) )
X		{
X		case PPM_TYPE:
X		uniform = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( uniform <= EPSILON )
X		    rayleigh = INFINITY;
X		else
X		    rayleigh = sqrt( -2.0 * log( uniform ) );
X		uniform = ( random( ) & RANDOM_MASK ) / arand; 
X		gauss = rayleigh * cos( 2.0 * M_PI * uniform );
X		rtemp = PPM_GETR( *xP ) + ( PPM_GETR( *xP ) * mgsigma * gauss );
X
X		uniform = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( uniform <= EPSILON )
X		    rayleigh = INFINITY;
X		else
X		    rayleigh = sqrt( -2.0 * log( uniform ) );
X		uniform = ( random( ) & RANDOM_MASK ) / arand; 
X		gauss = rayleigh * cos( 2.0 * M_PI * uniform );
X		gtemp = PPM_GETG( *xP ) + ( PPM_GETG( *xP ) * mgsigma * gauss );
X
X		uniform = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( uniform <= EPSILON )
X		    rayleigh = INFINITY;
X		else
X		    rayleigh = sqrt( -2.0 * log( uniform ) );
X		uniform = ( random( ) & RANDOM_MASK ) / arand; 
X		gauss = rayleigh * cos( 2.0 * M_PI * uniform );
X		btemp = PPM_GETB( *xP ) + ( PPM_GETB( *xP ) * mgsigma * gauss );
X
X		CHECK_RED;
X		CHECK_GREEN;
X		CHECK_BLUE;
X		PPM_ASSIGN( *xP, r, g, b );
X		break;
X
X		default:
X		uniform = ( random( ) & RANDOM_MASK ) / arand; 
X		if ( uniform <= EPSILON )
X		    rayleigh = INFINITY;
X		else
X		    rayleigh = sqrt( -2.0 * log( uniform ) );
X		uniform = ( random( ) & RANDOM_MASK ) / arand; 
X		gauss = rayleigh * cos( 2.0 * M_PI * uniform );
X		gtemp = PNM_GET1( *xP ) + ( PNM_GET1( *xP ) * mgsigma * gauss );
X
X		CHECK_GRAY;
X		PNM_ASSIGN1( *xP, g );
X		break;
X		}
X	    }
X	pnm_writepnmrow( stdout, xelrow, cols, maxval, format, forceplain );
X	}
X    pnm_freerow( xelrow );
X    }
X
X/* Poisson Noise
X**
X*/
X
Xstatic void
Xpoisson_noise( FILE *ifp, float lambda )
X    {
X    register xel *xP;
X    register int col;
X    xel *xelrow;
X    xelval r, g, b;
X    xelval maxval;
X    int rows, cols, format, row;
X    double rtemp, gtemp, btemp;
X    double x, x1;
X    float rr;
X    int k;
X
X    int forceplain = 0;
X
X    pnm_readpnminit( ifp, &cols, &rows, &maxval, &format );
X    xelrow = pnm_allocrow( cols );
X
X    pnm_writepnminit( stdout, cols, rows, maxval, format, forceplain );
X
X    for ( row = 0; row < rows; ++row )
X	{
X	pnm_readpnmrow( ifp, xelrow, cols, maxval, format );
X	for ( col = 0, xP = xelrow; col < cols; ++col, ++xP )
X	    {
X	    switch ( PNM_FORMAT_TYPE( format ) )
X		{
X		case PPM_TYPE:
X		x = lambda * PPM_GETR( *xP );
X		x1 = exp( -x );
X		rr = 1.0;
X		k = 0;
X		rr = rr * ( ( random( ) & RANDOM_MASK ) / arand );
X		while ( rr > x1 )
X		    {
X		    ++k;
X		    rr = rr * ( ( random( ) & RANDOM_MASK ) / arand );
X		    }
X		rtemp = k / lambda;
X
X		x = lambda * PPM_GETG( *xP );
X		x1 = exp( -x );
X		rr = 1.0;
X		k = 0;
X		rr = rr * ( ( random( ) & RANDOM_MASK ) / arand );
X		while ( rr > x1 )
X		    {
X		    ++k;
X		    rr = rr * ( ( random( ) & RANDOM_MASK ) / arand );
X		    }    
X		gtemp = k / lambda;
X
X		x = lambda * PPM_GETB( *xP );
X		x1 = exp( -x );
X		rr = 1.0;
X		k = 0;
X		rr = rr * ( ( random( ) & RANDOM_MASK ) / arand );
X		while ( rr > x1 )
X		    {
X		    ++k;
X		    rr = rr * ( ( random( ) & RANDOM_MASK ) / arand );
X		    }
X		btemp = k / lambda;
X
X		CHECK_RED;
X		CHECK_GREEN;
X		CHECK_BLUE;
X		PPM_ASSIGN( *xP, r, g, b );
X		break;
X
X		default:
X		x = lambda * PNM_GET1( *xP );
X		x1 = exp( -x );
X		rr = 1.0;
X		k = 0;
X		rr = rr * ( ( random( ) & RANDOM_MASK ) / arand );
X		while ( rr > x1 )
X		    {
X		    ++k;
X		    rr = rr * ( ( random( ) & RANDOM_MASK ) / arand );
X		    }
X		gtemp = k / lambda;
X
X		CHECK_GRAY;
X		PNM_ASSIGN1( *xP, g );
X		break;
X		}
X	    }
X	pnm_writepnmrow( stdout, xelrow, cols, maxval, format, forceplain );
X	}
X    pnm_freerow( xelrow );
X    }
SHAR_EOF
chmod 0644 pnmaddnoise.c || echo "restore of pnmaddnoise.c fails"
exit 0


