path: root/editor/pamaddnoise.c

/*
** 
** Add gaussian, multiplicative gaussian, impulse, laplacian or 
** poisson noise to a portable anymap.
** 
** Version 1.0  November 1995
**
** Copyright (C) 1995 by Mike Burns (burns@cac.psu.edu)
**
** Adapted to Netpbm 2005.08.09, by Bryan Henderson
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/

/* References
** ----------
** "Adaptive Image Restoration in Signal-Dependent Noise" by R. Kasturi
** Institute for Electronic Science, Texas Tech University  1982
**
** "Digital Image Processing Algorithms" by Ioannis Pitas
** Prentice Hall, 1993  ISBN 0-13-145814-0
*/

#define _XOPEN_SOURCE 500  /* get M_PI in math.h */

#include <math.h>

#include "pm_c_util.h"
#include "pam.h"

#define RANDOM_MASK 0x7FFF  /* only compare lower 15 bits.  Stupid PCs. */

static double const EPSILON = 1.0e-5;
static double const arand = 32767.0;      /* 2^15-1 in case stoopid computer */

enum noiseType {
    GAUSSIAN,
    IMPULSE,  /* aka salt and pepper noise */
    LAPLACIAN,
    MULTIPLICATIVE_GAUSSIAN,
    POISSON,
    MAX_NOISE_TYPES
};



static void
gaussian_noise(sample   const maxval,
               sample   const origSample,
               sample * const newSampleP,
               float    const sigma1,
               float    const sigma2) {
/*----------------------------------------------------------------------------
   Add Gaussian noise.

   Based on Kasturi/Algorithms of the ACM
-----------------------------------------------------------------------------*/

    double x1, x2, xn, yn;
    double rawNewSample;

    x1 = (rand() & RANDOM_MASK) / arand; 

    if (x1 == 0.0)
        x1 = 1.0;
    x2 = (rand() & RANDOM_MASK) / arand;
    xn = sqrt(-2.0 * log(x1)) * cos(2.0 * M_PI * x2);
    yn = sqrt(-2.0 * log(x1)) * sin(2.0 * M_PI * x2);
    
    rawNewSample =
        origSample + (sqrt((double) origSample) * sigma1 * xn) + (sigma2 * yn);

    *newSampleP = MAX(MIN((int)rawNewSample, maxval), 0);
}



static void
impulse_noise(sample   const maxval,
              sample   const origSample,
              sample * const newSampleP,
              float    const tolerance) {
/*----------------------------------------------------------------------------
   Add impulse (salt and pepper) noise
-----------------------------------------------------------------------------*/

    double const low_tol  = tolerance / 2.0;
    double const high_tol = 1.0 - (tolerance / 2.0);
    double const sap = (rand() & RANDOM_MASK) / arand; 

    if (sap < low_tol) 
        *newSampleP = 0;
    else if ( sap >= high_tol )
        *newSampleP = maxval;
}



static void
laplacian_noise(sample   const maxval,
                double   const infinity,
                sample   const origSample,
                sample * const newSampleP,
                float    const lsigma) {
/*----------------------------------------------------------------------------
   Add Laplacian noise

   From Pitas' book.
-----------------------------------------------------------------------------*/
    double const u = (rand() & RANDOM_MASK) / arand; 
                
    double rawNewSample;

    if (u <= 0.5) {
        if (u <= EPSILON)
            rawNewSample = origSample - infinity;
        else
            rawNewSample = origSample + lsigma * log(2.0 * u);
    } else {
        double const u1 = 1.0 - u;
        if (u1 <= 0.5 * EPSILON)
            rawNewSample = origSample + infinity;
        else
            rawNewSample = origSample - lsigma * log(2.0 * u1);
    }
    *newSampleP = MIN(MAX((int)rawNewSample, 0), maxval);
}



static void
multiplicative_gaussian_noise(sample   const maxval,
                              double   const infinity,
                              sample   const origSample,
                              sample * const newSampleP,
                              float    const mgsigma) {
/*----------------------------------------------------------------------------
   Add multiplicative Gaussian noise

   From Pitas' book.
-----------------------------------------------------------------------------*/
    double rayleigh, gauss;
    double rawNewSample;

    {
        double const uniform = (rand() & RANDOM_MASK) / arand; 
        if (uniform <= EPSILON)
            rayleigh = infinity;
        else
            rayleigh = sqrt(-2.0 * log( uniform));
    }
    {
        double const uniform = (rand() & RANDOM_MASK) / arand; 
        gauss = rayleigh * cos(2.0 * M_PI * uniform);
    }
    rawNewSample = origSample + (origSample * mgsigma * gauss);

    *newSampleP = MIN(MAX((int)rawNewSample, 0), maxval);
}



static void
poisson_noise(sample   const maxval,
              sample   const origSample,
              sample * const newSampleP,
              float    const lambda) {
/*----------------------------------------------------------------------------
   Add Poisson noise
-----------------------------------------------------------------------------*/
    double const x  = lambda * origSample;
    double const x1 = exp(-x);

    double rawNewSample;
    float rr;
    unsigned int k;

    rr = 1.0;  /* initial value */
    k = 0;     /* initial value */
    rr = rr * ((rand() & RANDOM_MASK) / arand);
    while (rr > x1) {
        ++k;
        rr = rr * ((rand() & RANDOM_MASK) / arand);
    }
    rawNewSample = k / lambda;

    *newSampleP = MIN(MAX((int)rawNewSample, 0), maxval);
}



int 
main(int argc, char * argv[]) {

    FILE * ifP;
    struct pam inpam;
    struct pam outpam;
    tuple * tuplerow;
    const tuple * newtuplerow;
    unsigned int row;
    double infinity;

    int argn;
    const char * inputFilename;
    int noise_type;
    unsigned int seed;
    int i;
    const char * const usage = "[-type noise_type] [-lsigma x] [-mgsigma x] "
        "[-sigma1 x] [-sigma2 x] [-lambda x] [-seed n] "
        "[-tolerance ratio] [pgmfile]";

    const char * const noise_name[] = { 
        "gaussian",
        "impulse",
        "laplacian",
        "multiplicative_gaussian",
        "poisson"
    };
    int const noise_id[] = { 
        GAUSSIAN,
        IMPULSE,
        LAPLACIAN,
        MULTIPLICATIVE_GAUSSIAN,
        POISSON
    };
    /* minimum number of characters to match noise name for pm_keymatch() */
    int const noise_compare[] = {
        1,
        1,
        1,
        1,
        1
    };

    /* define default values for configurable options */
    float lambda = 0.05;        
    float lsigma = 10.0;
    float mgsigma = 0.5;
    float sigma1 = 4.0;
    float sigma2 = 20.0;
    float tolerance = 0.10;

    pnm_init(&argc, argv);

    seed = pm_randseed();
    noise_type = GAUSSIAN;

    argn = 1;
    while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
    {
        if ( pm_keymatch( argv[argn], "-lambda", 3 ) )
        {
            ++argn;
            if ( argn >= argc )
            {
                pm_message( 
                    "incorrect number of arguments for -lambda option" );
                pm_usage( usage );
            }
            else if ( argv[argn][0] == '-' )
            {
                pm_message( "invalid argument to -lambda option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            lambda = atof( argv[argn] );
        }
        else if ( pm_keymatch( argv[argn], "-lsigma", 3 ) )
        {
            ++argn;
            if ( argn >= argc )
            {
                pm_message( 
                    "incorrect number of arguments for -lsigma option" );
                pm_usage( usage );
            }
            else if ( argv[argn][0] == '-' )
            {
                pm_message( "invalid argument to -lsigma option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            lsigma = atof( argv[argn] );
        }
        else if ( pm_keymatch( argv[argn], "-mgsigma", 2 ) )
        {
            ++argn;
            if ( argn >= argc )
            {
                pm_message( 
                    "incorrect number of arguments for -mgsigma option" );
                pm_usage( usage );
            }
            else if ( argv[argn][0] == '-' )
            {
                pm_message( "invalid argument to -mgsigma option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            mgsigma = atof( argv[argn] );
        }
        else if ( pm_keymatch( argv[argn], "-seed", 3 ) )
        {
            ++argn;
            if ( argn >= argc )
            {
                pm_message( "incorrect number of arguments for -seed option" );
                pm_usage( usage );
            }
            else if ( argv[argn][0] == '-' )
            {
                pm_message( "invalid argument to -seed option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            seed = atoi(argv[argn]);
        }
        else if ( pm_keymatch( argv[argn], "-sigma1", 7 ) ||
                  pm_keymatch( argv[argn], "-s1", 3 ) )
        {
            ++argn;
            if ( argn >= argc )
            {
                pm_message( 
                    "incorrect number of arguments for -sigma1 option" );
                pm_usage( usage );
            }
            else if ( argv[argn][0] == '-' )
            {
                pm_message( "invalid argument to -sigma1 option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            sigma1 = atof( argv[argn] );
        }
        else if ( pm_keymatch( argv[argn], "-sigma2", 7 ) ||
                  pm_keymatch( argv[argn], "-s2", 3 ) )
        {
            ++argn;
            if ( argn >= argc )
            {
                pm_message( 
                    "incorrect number of arguments for -sigma2 option" );
                pm_usage( usage );
            }
            else if ( argv[argn][0] == '-' )
            {
                pm_message( "invalid argument to -sigma2 option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            sigma2 = atof( argv[argn] );
        }
        else if ( pm_keymatch( argv[argn], "-tolerance", 3 ) )
        {
            ++argn;
            if ( argn >= argc )
            {
                pm_message( 
                    "incorrect number of arguments for -tolerance option" );
                pm_usage( usage );
            }
            else if ( argv[argn][0] == '-' )
            {
                pm_message( "invalid argument to -tolerance option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            tolerance = atof( argv[argn] );
        }
        else if ( pm_keymatch( argv[argn], "-type", 3 ) )
        {
            ++argn;
            if ( argn >= argc )
            {
                pm_message( "incorrect number of arguments for -type option" );
                pm_usage( usage );
            }
            else if ( argv[argn][0] == '-' )
            {
                pm_message( "invalid argument to -type option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            /* search through list of valid noise types and compare */
            i = 0;
            while ( ( i < MAX_NOISE_TYPES ) && 
                    !pm_keymatch( argv[argn], 
                                  noise_name[i], noise_compare[i] ) )
                ++i;
            if ( i >= MAX_NOISE_TYPES )
            {
                pm_message( "invalid argument to -type option: %s", 
                            argv[argn] );
                pm_usage( usage );
            }
            noise_type = noise_id[i];
        }
        else
            pm_usage( usage );
        ++argn;
    }

    if ( argn < argc )
    {
        inputFilename = argv[argn];
        argn++;
    }
    else
        inputFilename = "-";

    if ( argn != argc )
        pm_usage( usage );

    srand(seed);

    ifP = pm_openr(inputFilename);

    pnm_readpaminit(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type));

    outpam = inpam;
    outpam.file = stdout;

    pnm_writepaminit(&outpam);

    tuplerow = pnm_allocpamrow(&inpam);
    newtuplerow = pnm_allocpamrow(&inpam);
    infinity = (double) inpam.maxval;
    
    for (row = 0; row < inpam.height; ++row) {
        unsigned int col;
        pnm_readpamrow(&inpam, tuplerow);
        for (col = 0; col < inpam.width; ++col) {
            unsigned int plane;
            for (plane = 0; plane < inpam.depth; ++plane) {
                switch (noise_type) {
                case GAUSSIAN:
                    gaussian_noise(inpam.maxval,
                                   tuplerow[col][plane],
                                   &newtuplerow[col][plane],
                                   sigma1, sigma2);
                    break;
                    
                case IMPULSE:
                    impulse_noise(inpam.maxval,
                                  tuplerow[col][plane],
                                  &newtuplerow[col][plane],
                                  tolerance);
                   break;
                    
                case LAPLACIAN:
                    laplacian_noise(inpam.maxval, infinity,
                                    tuplerow[col][plane],
                                    &newtuplerow[col][plane],
                                    lsigma);
                    break;
                    
                case MULTIPLICATIVE_GAUSSIAN:
                    multiplicative_gaussian_noise(inpam.maxval, infinity,
                                                  tuplerow[col][plane],
                                                  &newtuplerow[col][plane],
                                                  mgsigma);
                    break;
                    
                case POISSON:
                    poisson_noise(inpam.maxval,
                                  tuplerow[col][plane],
                                  &newtuplerow[col][plane],
                                  lambda);
                    break;

                }
            }
        }
        pnm_writepamrow(&outpam, newtuplerow);
    }
    pnm_freepamrow(newtuplerow);
    pnm_freepamrow(tuplerow);

    return 0;
}