diff options
Diffstat (limited to 'editor/pamaddnoise.c')
| -rw-r--r-- | editor/pamaddnoise.c | 236 |
1 files changed, 138 insertions, 98 deletions
diff --git a/editor/pamaddnoise.c b/editor/pamaddnoise.c index 20c68d99..b662aa96 100644 --- a/editor/pamaddnoise.c +++ b/editor/pamaddnoise.c @@ -33,20 +33,20 @@ #include "pm_c_util.h" #include "mallocvar.h" +#include "rand.h" #include "shhopt.h" #include "pm_gamma.h" #include "pam.h" static double const EPSILON = 1.0e-5; - - -static double -rand1() { - - return (double)rand()/RAND_MAX; -} - +static double const SIGMA1_DEFAULT = 4.0; +static double const SIGMA2_DEFAULT = 20.0; +static double const MGSIGMA_DEFAULT = 0.5; +static double const LSIGMA_DEFAULT = 10.0; +static double const TOLERANCE_DEFAULT = 0.10; +static double const SALT_RATIO_DEFAULT = 0.5; +static double const LAMBDA_DEFAULT = 12.0; enum NoiseType { @@ -67,6 +67,7 @@ struct CmdlineInfo { enum NoiseType noiseType; + unsigned int seedSpec; unsigned int seed; float lambda; @@ -75,6 +76,7 @@ struct CmdlineInfo { float sigma1; float sigma2; float tolerance; + float saltRatio; }; @@ -119,8 +121,8 @@ parseCommandLine(int argc, const char ** const argv, unsigned int option_def_index; - unsigned int typeSpec, seedSpec, lambdaSpec, lsigmaSpec, mgsigmaSpec, - sigma1Spec, sigma2Spec, toleranceSpec; + unsigned int typeSpec, lambdaSpec, lsigmaSpec, mgsigmaSpec, + sigma1Spec, sigma2Spec, toleranceSpec, saltRatioSpec; const char * type; @@ -128,21 +130,23 @@ parseCommandLine(int argc, const char ** const argv, option_def_index = 0; /* incremented by OPTENT3 */ OPTENT3(0, "type", OPT_STRING, &type, - &typeSpec, 0); + &typeSpec, 0); OPTENT3(0, "seed", OPT_UINT, &cmdlineP->seed, - &seedSpec, 0); + &cmdlineP->seedSpec, 0); OPTENT3(0, "lambda", OPT_FLOAT, &cmdlineP->lambda, - &lambdaSpec, 0); + &lambdaSpec, 0); OPTENT3(0, "lsigma", OPT_FLOAT, &cmdlineP->lsigma, - &lsigmaSpec, 0); + &lsigmaSpec, 0); OPTENT3(0, "mgsigma", OPT_FLOAT, &cmdlineP->mgsigma, - &mgsigmaSpec, 0); + &mgsigmaSpec, 0); OPTENT3(0, "sigma1", OPT_FLOAT, &cmdlineP->sigma1, - &sigma1Spec, 0); + &sigma1Spec, 0); OPTENT3(0, "sigma2", OPT_FLOAT, &cmdlineP->sigma2, - &sigma2Spec, 0); + &sigma2Spec, 0); OPTENT3(0, "tolerance", OPT_FLOAT, &cmdlineP->tolerance, - &toleranceSpec, 0); + &toleranceSpec, 0); + OPTENT3(0, "salt", OPT_FLOAT, &cmdlineP->saltRatio, + &saltRatioSpec, 0); opt.opt_table = option_def; opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ @@ -156,44 +160,84 @@ parseCommandLine(int argc, const char ** const argv, else cmdlineP->noiseType = typeFmName(type); - if (sigma1Spec && cmdlineP->noiseType != NOISETYPE_GAUSSIAN) - pm_error("-sigma1 is valid only with -type=gaussian"); + if (sigma1Spec) { + if (cmdlineP->noiseType != NOISETYPE_GAUSSIAN) + pm_error("-sigma1 is valid only with -type=gaussian"); + else if (cmdlineP->sigma1 < 0) + pm_error("-sigma1 value must be non-negative. You specified %f", + cmdlineP->sigma1); + } + + if (sigma2Spec) { + if (cmdlineP->noiseType != NOISETYPE_GAUSSIAN) + pm_error("-sigma2 is valid only with -type=gaussian"); + else if (cmdlineP->sigma2 < 0) + pm_error("-sigma2 value must be non-negative. You specified %f", + cmdlineP->sigma2); + } - if (sigma2Spec && cmdlineP->noiseType != NOISETYPE_GAUSSIAN) - pm_error("-sigma2 is valid only with -type=gaussian"); + if (mgsigmaSpec) { + if (cmdlineP->noiseType != NOISETYPE_MULTIPLICATIVE_GAUSSIAN) + pm_error("-mgsigma is valid only with -type=multiplicative_guassian"); + else if (cmdlineP->mgsigma < 0) + pm_error("-mgsigma value must be non-negative. You specified %f", + cmdlineP->mgsigma); + } - if (mgsigmaSpec && - cmdlineP->noiseType != NOISETYPE_MULTIPLICATIVE_GAUSSIAN) - pm_error("-mgsigma is valid only with -type=multiplicative_guassian"); + if (toleranceSpec) { + if (cmdlineP->noiseType != NOISETYPE_IMPULSE) + pm_error("-tolerance is valid only with -type=impulse"); + else if (cmdlineP->tolerance < 0 || cmdlineP->tolerance > 1.0) + pm_error("-tolerance value must be between 0.0 and 1.0. " + "You specified %f", cmdlineP->tolerance); + } - if (toleranceSpec && cmdlineP->noiseType != NOISETYPE_IMPULSE) - pm_error("-tolerance is valid only with -type=impulse"); + if (saltRatioSpec) { + if (cmdlineP->noiseType != NOISETYPE_IMPULSE) + pm_error("-salt is valid only with -type=impulse"); + else if (cmdlineP->saltRatio < 0 || cmdlineP->saltRatio > 1.0) + pm_error("-salt value must be between 0.0 and 1.0. " + "You specified %f", cmdlineP->saltRatio); + } - if (lsigmaSpec && cmdlineP->noiseType != NOISETYPE_LAPLACIAN) + if (lsigmaSpec) { + if (cmdlineP->noiseType != NOISETYPE_LAPLACIAN) pm_error("-lsigma is valid only with -type=laplacian"); + else if (cmdlineP->lsigma <= 0) + pm_error("-lsigma value must be positive. You specified %f", + cmdlineP->lsigma); + } - if (lambdaSpec && cmdlineP->noiseType != NOISETYPE_POISSON) + if (lambdaSpec) { + if (cmdlineP->noiseType != NOISETYPE_POISSON) pm_error("-lambda is valid only with -type=poisson"); + else if (cmdlineP->lambda <= 0) + pm_error("-lambda value must be positive. You specified %f", + cmdlineP->lambda); + } if (!lambdaSpec) - cmdlineP->lambda = 12.0; + cmdlineP->lambda = LAMBDA_DEFAULT; if (!lsigmaSpec) - cmdlineP->lsigma = 10.0; + cmdlineP->lsigma = LSIGMA_DEFAULT; if (!mgsigmaSpec) - cmdlineP->mgsigma = 0.5; + cmdlineP->mgsigma = MGSIGMA_DEFAULT; if (!sigma1Spec) - cmdlineP->sigma1 = 4.0; + cmdlineP->sigma1 = SIGMA1_DEFAULT; if (!sigma2Spec) - cmdlineP->sigma2 = 20.0; + cmdlineP->sigma2 = SIGMA2_DEFAULT; if (!toleranceSpec) - cmdlineP->tolerance = 0.10; + cmdlineP->tolerance = TOLERANCE_DEFAULT; + + if (!saltRatioSpec) + cmdlineP->saltRatio = SALT_RATIO_DEFAULT; - if (!seedSpec) + if (!cmdlineP->seedSpec) cmdlineP->seed = pm_randseed(); if (argc-1 > 1) @@ -211,30 +255,25 @@ parseCommandLine(int argc, const char ** const argv, static void -addGaussianNoise(sample const maxval, - sample const origSample, - sample * const newSampleP, - float const sigma1, - float const sigma2) { +addGaussianNoise(sample const maxval, + sample const origSample, + sample * const newSampleP, + float const sigma1, + float const sigma2, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- Add Gaussian noise. Based on Kasturi/Algorithms of the ACM -----------------------------------------------------------------------------*/ - double x1, x2, xn, yn; + double grnd1, grnd2; /* Gaussian random numbers. mean=0 sigma=1 */ double rawNewSample; - x1 = rand1(); - - if (x1 == 0.0) - x1 = 1.0; - x2 = rand1(); - xn = sqrt(-2.0 * log(x1)) * cos(2.0 * M_PI * x2); - yn = sqrt(-2.0 * log(x1)) * sin(2.0 * M_PI * x2); + pm_gaussrand2(randStP, &grnd1, &grnd2); rawNewSample = - origSample + (sqrt((double) origSample) * sigma1 * xn) + (sigma2 * yn); + origSample + (sqrt((double) origSample) * sigma1 * grnd1) + (sigma2 * grnd2); *newSampleP = MAX(MIN((int)rawNewSample, maxval), 0); } @@ -242,40 +281,42 @@ addGaussianNoise(sample const maxval, static void -addImpulseNoise(sample const maxval, - sample const origSample, - sample * const newSampleP, - float const tolerance) { +addImpulseNoise(sample const maxval, + sample const origSample, + sample * const newSampleP, + float const tolerance, + double const saltRatio, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- 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 = rand1(); + double const pepperRatio = 1.0 - saltRatio; + double const loTolerance = tolerance * pepperRatio; + double const hiTolerance = 1.0 - tolerance * saltRatio; + double const sap = pm_drand(randStP); - if (sap < low_tol) - *newSampleP = 0; - else if ( sap >= high_tol ) - *newSampleP = maxval; - else - *newSampleP = origSample; + *newSampleP = + sap < loTolerance ? 0 : + sap >= hiTolerance? maxval : + origSample; } static void -addLaplacianNoise(sample const maxval, - double const infinity, - sample const origSample, - sample * const newSampleP, - float const lsigma) { +addLaplacianNoise(sample const maxval, + double const infinity, + sample const origSample, + sample * const newSampleP, + float const lsigma, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- Add Laplacian noise From Pitas' book. -----------------------------------------------------------------------------*/ - double const u = rand1(); + double const u = pm_drand(randStP); double rawNewSample; @@ -297,31 +338,21 @@ addLaplacianNoise(sample const maxval, static void -addMultiplicativeGaussianNoise(sample const maxval, - double const infinity, - sample const origSample, - sample * const newSampleP, - float const mgsigma) { +addMultiplicativeGaussianNoise(sample const maxval, + double const infinity, + sample const origSample, + sample * const newSampleP, + float const mgsigma, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- Add multiplicative Gaussian noise From Pitas' book. -----------------------------------------------------------------------------*/ - double rayleigh, gauss; + double rawNewSample; - { - double const uniform = rand1(); - if (uniform <= EPSILON) - rayleigh = infinity; - else - rayleigh = sqrt(-2.0 * log( uniform)); - } - { - double const uniform = rand1(); - gauss = rayleigh * cos(2.0 * M_PI * uniform); - } - rawNewSample = origSample + (origSample * mgsigma * gauss); + rawNewSample = origSample + (origSample * mgsigma * pm_gaussrand(randStP)); *newSampleP = MIN(MAX((int)rawNewSample, 0), maxval); } @@ -363,10 +394,11 @@ poissonPmf(double const lambda, static void -addPoissonNoise(struct pam * const pamP, - sample const origSample, - sample * const newSampleP, - float const lambdaOfMaxval) { +addPoissonNoise(struct pam * const pamP, + sample const origSample, + sample * const newSampleP, + float const lambdaOfMaxval, + struct pm_randSt * const randStP) { /*---------------------------------------------------------------------------- Add Poisson noise -----------------------------------------------------------------------------*/ @@ -376,7 +408,7 @@ addPoissonNoise(struct pam * const pamP, double const lambda = origSampleIntensity * lambdaOfMaxval; - double const u = rand1(); + double const u = pm_drand(randStP); /* We now apply the inverse CDF (cumulative distribution function) of the Poisson distribution to uniform random variable 'u' to get a Poisson @@ -416,12 +448,14 @@ main(int argc, const char ** argv) { const tuple * newtuplerow; unsigned int row; double infinity; + struct pm_randSt randSt; pm_proginit(&argc, argv); parseCommandLine(argc, argv, &cmdline); - srand(cmdline.seed); + pm_randinit(&randSt); + pm_srand2(&randSt, cmdline.seedSpec, cmdline.seed); ifP = pm_openr(cmdline.inputFileName); @@ -448,35 +482,40 @@ main(int argc, const char ** argv) { addGaussianNoise(inpam.maxval, tuplerow[col][plane], &newtuplerow[col][plane], - cmdline.sigma1, cmdline.sigma2); + cmdline.sigma1, cmdline.sigma2, + &randSt); break; case NOISETYPE_IMPULSE: addImpulseNoise(inpam.maxval, tuplerow[col][plane], &newtuplerow[col][plane], - cmdline.tolerance); + cmdline.tolerance, cmdline.saltRatio, + &randSt); break; case NOISETYPE_LAPLACIAN: addLaplacianNoise(inpam.maxval, infinity, tuplerow[col][plane], &newtuplerow[col][plane], - cmdline.lsigma); + cmdline.lsigma, + &randSt); break; case NOISETYPE_MULTIPLICATIVE_GAUSSIAN: addMultiplicativeGaussianNoise(inpam.maxval, infinity, tuplerow[col][plane], &newtuplerow[col][plane], - cmdline.mgsigma); + cmdline.mgsigma, + &randSt); break; case NOISETYPE_POISSON: addPoissonNoise(&inpam, tuplerow[col][plane], &newtuplerow[col][plane], - cmdline.lambda); + cmdline.lambda, + &randSt); break; } @@ -484,6 +523,7 @@ main(int argc, const char ** argv) { } pnm_writepamrow(&outpam, newtuplerow); } + pm_randterm(&randSt); pnm_freepamrow(newtuplerow); pnm_freepamrow(tuplerow); |