1## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
2## Copyright (C) 2004 Stefan van der Walt <stefan@sun.ac.za>
3## Copyright (C) 2012 Carlo de Falco
4## Copyright (C) 2012 Carnë Draug <carandraug@octave.org>
5##
6## This program is free software; you can redistribute it and/or modify it under
7## the terms of the GNU General Public License as published by the Free Software
8## Foundation; either version 3 of the License, or (at your option) any later
9## version.
10##
11## This program is distributed in the hope that it will be useful, but WITHOUT
12## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14## details.
15##
16## You should have received a copy of the GNU General Public License along with
17## this program; if not, see <http://www.gnu.org/licenses/>.
18
19## -*- texinfo -*-
20## @deftypefn  {Function File} {} imnoise (@var{A}, @var{type})
21## @deftypefnx {Function File} {} imnoise (@dots{}, @var{options})
22## Add noise to image.
23##
24## @end deftypefn
25## @deftypefn {Function File} {} imnoise (@var{A}, "gaussian", @var{mean}, @var{variance})
26## Additive gaussian noise with @var{mean} and @var{variance} defaulting to 0
27## and 0.01.
28##
29## @end deftypefn
30## @deftypefn {Function File} {} imnoise (@var{A}, "poisson")
31## Creates poisson noise in the image using the intensity value of each pixel as
32## mean.
33##
34## @end deftypefn
35## @deftypefn {Function File} {} imnoise (@var{A}, "salt & pepper", @var{density})
36## Create "salt and pepper"/"lost pixels" in @var{density}*100 percent of the
37## image.  @var{density} defaults to 0.05.
38##
39## @end deftypefn
40## @deftypefn {Function File} {} imnoise (@var{A}, "speckle", @var{variance})
41## Multiplicative gaussian noise with @var{B} = @var{A} + @var{A} * noise with
42## mean 0 and @var{variance} defaulting to 0.04.
43##
44## @seealso{rand, randn, randp}
45## @end deftypefn
46
47function A = imnoise (A, stype, a, b)
48  ## we do not set defaults right at the start because they are different
49  ## depending on the method used to generate noise
50
51  if (nargin < 2 || nargin > 4)
52    print_usage;
53  elseif (! isimage (A))
54    error ("imnoise: first argument must be an image.");
55  elseif (! ischar (stype))
56    error ("imnoise: second argument must be a string with name of noise type.");
57  endif
58
59  in_class  = class (A);
60  fix_class = false;      # for cases when we need to use im2double
61
62  switch (lower (stype))
63    case "poisson"
64      switch (in_class)
65        case ("double")
66          A = randp (A * 1e12) / 1e12;
67        case ("single")
68          A = single (randp (A * 1e6) / 1e6);
69        case {"uint8", "uint16"}
70          A = cast (randp (A), in_class);
71        otherwise
72          A = imnoise (im2double (A), "poisson");
73          fix_class = true;
74      endswitch
75
76    case "gaussian"
77      A         = im2double (A);
78      fix_class = true;
79      if (nargin < 3), a = 0.00; endif
80      if (nargin < 4), b = 0.01; endif
81      A = A + (a + randn (size (A)) * sqrt (b));
82      ## Variance of Gaussian data with mean 0 is E[X^2]
83
84    case {"salt & pepper", "salt and pepper"}
85      if (nargin < 3), a = 0.05; endif
86      noise = rand (size (A));
87      if (isfloat (A))
88        black = 0;
89        white = 1;
90      else
91        black = intmin (in_class);
92        white = intmax (in_class);
93      endif
94      A(noise <= a/2)   = black;
95      A(noise >= 1-a/2) = white;
96
97    case "speckle"
98      A         = im2double (A);
99      fix_class = true;
100      if (nargin < 3), a = 0.04; endif
101      A = A .* (1 + randn (size (A)) * sqrt (a));
102
103    otherwise
104      error ("imnoise: unknown or unimplemented type of noise `%s'", stype);
105  endswitch
106
107  if (fix_class)
108    A = imcast (A, in_class);
109  elseif (isfloat (A))
110    ## this includes not even cases where the noise made it go outside of the
111    ## [0 1] range, but also images that were already originally outside that
112    ## range. This is by design and matlab compatibility. And we do this after
113    ## fixing class because the imcast functions already take care of such
114    ## adjustment
115    A(A < 0) = cast (0, class (A));
116    A(A > 1) = cast (1, class (A));
117  endif
118
119endfunction
120
121%!assert(var(imnoise(ones(10)/2,'gaussian')(:)),0.01,0.005) # probabilistic
122%!assert(length(find(imnoise(ones(10)/2,'salt & pepper')~=0.5)),5,10) # probabilistic
123%!assert(var(imnoise(ones(10)/2,'speckle')(:)),0.01,0.005) # probabilistic
124
125%!test
126%! A = imnoise (.5 * ones (100), 'poisson');
127%! assert (class (A), 'double')
128%!test
129%! A = imnoise (.5 * ones (100, 'single'), 'poisson');
130%! assert (class (A), 'single')
131%!test
132%! A = imnoise (128 * ones (100, 'uint8'), 'poisson');
133%! assert (class (A), 'uint8')
134%!test
135%! A = imnoise (256 * ones (100, 'uint16'), 'poisson');
136%! assert (class (A), 'uint16')
137
138%!demo
139%!  A = imnoise (2^7 * ones (100, 'uint8'), 'poisson');
140%!  subplot (2, 2, 1)
141%!  imshow (A)
142%!  title ('uint8 image with poisson noise')
143%!  A = imnoise (2^15 * ones (100, 'uint16'), 'poisson');
144%!  subplot (2, 2, 2)
145%!  imshow (A)
146%!  title ('uint16 image with poisson noise')
147%!  A = imnoise (.5 * ones (100), 'poisson');
148%!  subplot (2, 2, 3)
149%!  imshow (A)
150%!  title ('double image with poisson noise')
151%!  A = imnoise (.5 * ones (100, 'single'), 'poisson');
152%!  subplot (2, 2, 4)
153%!  imshow (A)
154%!  title ('single image with poisson noise')
155