MagickCore  7.1.1-43
Convert, Edit, Or Compose Bitmap Images
gem.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % GGGG EEEEE M M %
7 % G E MM MM %
8 % G GG EEE M M M %
9 % G G E M M %
10 % GGGG EEEEE M M %
11 % %
12 % %
13 % Graphic Gems - Graphic Support Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % August 1996 %
18 % %
19 % %
20 % Copyright @ 1999 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39 
40 /*
41  Include declarations.
42 */
43 #include "MagickCore/studio.h"
44 #include "MagickCore/color-private.h"
45 #include "MagickCore/draw.h"
46 #include "MagickCore/gem.h"
47 #include "MagickCore/gem-private.h"
48 #include "MagickCore/image.h"
49 #include "MagickCore/image-private.h"
50 #include "MagickCore/log.h"
51 #include "MagickCore/memory_.h"
52 #include "MagickCore/pixel-accessor.h"
53 #include "MagickCore/quantum.h"
54 #include "MagickCore/quantum-private.h"
55 #include "MagickCore/random_.h"
56 #include "MagickCore/resize.h"
57 #include "MagickCore/transform.h"
58 #include "MagickCore/signature-private.h"
59 
60 /*
61 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62 % %
63 % %
64 % %
65 % E x p a n d A f f i n e %
66 % %
67 % %
68 % %
69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70 %
71 % ExpandAffine() computes the affine's expansion factor, i.e. the square root
72 % of the factor by which the affine transform affects area. In an affine
73 % transform composed of scaling, rotation, shearing, and translation, returns
74 % the amount of scaling.
75 %
76 % The format of the ExpandAffine method is:
77 %
78 % double ExpandAffine(const AffineMatrix *affine)
79 %
80 % A description of each parameter follows:
81 %
82 % o expansion: ExpandAffine returns the affine's expansion factor.
83 %
84 % o affine: A pointer the affine transform of type AffineMatrix.
85 %
86 */
87 MagickExport double ExpandAffine(const AffineMatrix *affine)
88 {
89  assert(affine != (const AffineMatrix *) NULL);
90  return(sqrt(fabs(affine->sx*affine->sy-affine->rx*affine->ry)));
91 }
92 
93 /*
94 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95 % %
96 % %
97 % %
98 % G e n e r a t e D i f f e r e n t i a l N o i s e %
99 % %
100 % %
101 % %
102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
103 %
104 % GenerateDifferentialNoise() generates differential noise.
105 %
106 % The format of the GenerateDifferentialNoise method is:
107 %
108 % double GenerateDifferentialNoise(RandomInfo *random_info,
109 % const Quantum pixel,const NoiseType noise_type,const double attenuate)
110 %
111 % A description of each parameter follows:
112 %
113 % o random_info: the random info.
114 %
115 % o pixel: noise is relative to this pixel value.
116 %
117 % o noise_type: the type of noise.
118 %
119 % o attenuate: attenuate the noise.
120 %
121 */
122 MagickPrivate double GenerateDifferentialNoise(RandomInfo *random_info,
123  const Quantum pixel,const NoiseType noise_type,const double attenuate)
124 {
125 #define SigmaUniform (attenuate*0.015625)
126 #define SigmaGaussian (attenuate*0.015625)
127 #define SigmaImpulse (attenuate*0.1)
128 #define SigmaLaplacian (attenuate*0.0390625)
129 #define SigmaMultiplicativeGaussian (attenuate*0.5)
130 #define SigmaPoisson (attenuate*12.5)
131 #define SigmaRandom (attenuate)
132 #define TauGaussian (attenuate*0.078125)
133 
134  double
135  alpha,
136  beta,
137  noise,
138  sigma;
139 
140  alpha=GetPseudoRandomValue(random_info);
141  switch (noise_type)
142  {
143  case UniformNoise:
144  default:
145  {
146  noise=(double) pixel+(double) QuantumRange*SigmaUniform*(alpha-0.5);
147  break;
148  }
149  case GaussianNoise:
150  {
151  double
152  gamma,
153  tau;
154 
155  if (fabs(alpha) < MagickEpsilon)
156  alpha=1.0;
157  beta=GetPseudoRandomValue(random_info);
158  gamma=sqrt(-2.0*log(alpha));
159  sigma=gamma*cos((double) (2.0*MagickPI*beta));
160  tau=gamma*sin((double) (2.0*MagickPI*beta));
161  noise=(double) pixel+sqrt((double) pixel)*SigmaGaussian*sigma+
162  (double) QuantumRange*TauGaussian*tau;
163  break;
164  }
165  case ImpulseNoise:
166  {
167  if (alpha < (SigmaImpulse/2.0))
168  noise=0.0;
169  else
170  if (alpha >= (1.0-(SigmaImpulse/2.0)))
171  noise=(double) QuantumRange;
172  else
173  noise=(double) pixel;
174  break;
175  }
176  case LaplacianNoise:
177  {
178  if (alpha <= 0.5)
179  {
180  if (alpha <= MagickEpsilon)
181  noise=(double) (pixel-QuantumRange);
182  else
183  noise=(double) pixel+(double) QuantumRange*SigmaLaplacian*
184  log(2.0*alpha)+0.5;
185  break;
186  }
187  beta=1.0-alpha;
188  if (beta <= (0.5*MagickEpsilon))
189  noise=(double) (pixel+QuantumRange);
190  else
191  noise=(double) pixel-(double) QuantumRange*SigmaLaplacian*
192  log(2.0*beta)+0.5;
193  break;
194  }
195  case MultiplicativeGaussianNoise:
196  {
197  sigma=1.0;
198  if (alpha > MagickEpsilon)
199  sigma=sqrt(-2.0*log(alpha));
200  beta=GetPseudoRandomValue(random_info);
201  noise=(double) pixel+(double) pixel*SigmaMultiplicativeGaussian*sigma*
202  cos((double) (2.0*MagickPI*beta))/2.0;
203  break;
204  }
205  case PoissonNoise:
206  {
207  double
208  poisson;
209 
210  ssize_t
211  i;
212 
213  poisson=exp(-SigmaPoisson*QuantumScale*(double) pixel);
214  for (i=0; alpha > poisson; i++)
215  {
216  beta=GetPseudoRandomValue(random_info);
217  alpha*=beta;
218  }
219  noise=(double) QuantumRange*i*PerceptibleReciprocal(SigmaPoisson);
220  break;
221  }
222  case RandomNoise:
223  {
224  noise=(double) QuantumRange*SigmaRandom*alpha;
225  break;
226  }
227  }
228  return(noise);
229 }
230 
231 /*
232 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
233 % %
234 % %
235 % %
236 % G e t O p t i m a l K e r n e l W i d t h %
237 % %
238 % %
239 % %
240 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
241 %
242 % GetOptimalKernelWidth() computes the optimal kernel radius for a convolution
243 % filter. Start with the minimum value of 3 pixels and walk out until we drop
244 % below the threshold of one pixel numerical accuracy.
245 %
246 % The format of the GetOptimalKernelWidth method is:
247 %
248 % size_t GetOptimalKernelWidth(const double radius,
249 % const double sigma)
250 %
251 % A description of each parameter follows:
252 %
253 % o width: GetOptimalKernelWidth returns the optimal width of a
254 % convolution kernel.
255 %
256 % o radius: the radius of the Gaussian, in pixels, not counting the center
257 % pixel.
258 %
259 % o sigma: the standard deviation of the Gaussian, in pixels.
260 %
261 */
262 MagickPrivate size_t GetOptimalKernelWidth1D(const double radius,
263  const double sigma)
264 {
265  double
266  alpha,
267  beta,
268  gamma,
269  normalize,
270  value;
271 
272  size_t
273  width;
274 
275  ssize_t
276  i,
277  j;
278 
279  if (IsEventLogging() != MagickFalse)
280  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
281  if (radius > MagickEpsilon)
282  return((size_t) (2.0*ceil(radius)+1.0));
283  gamma=fabs(sigma);
284  if (gamma <= MagickEpsilon)
285  return(3UL);
286  alpha=PerceptibleReciprocal(2.0*gamma*gamma);
287  beta=(double) PerceptibleReciprocal((double) MagickSQ2PI*gamma);
288  for (width=5; ; )
289  {
290  normalize=0.0;
291  j=(ssize_t) (width-1)/2;
292  for (i=(-j); i <= j; i++)
293  normalize+=exp(-((double) (i*i))*alpha)*beta;
294  value=exp(-((double) (j*j))*alpha)*beta/normalize;
295  if ((value < QuantumScale) || (value < MagickEpsilon))
296  break;
297  width+=2;
298  }
299  return((size_t) (width-2));
300 }
301 
302 MagickPrivate size_t GetOptimalKernelWidth2D(const double radius,
303  const double sigma)
304 {
305  double
306  alpha,
307  beta,
308  gamma,
309  normalize,
310  value;
311 
312  size_t
313  width;
314 
315  ssize_t
316  j,
317  u,
318  v;
319 
320  if (IsEventLogging() != MagickFalse)
321  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
322  if (radius > MagickEpsilon)
323  return((size_t) (2.0*ceil(radius)+1.0));
324  gamma=fabs(sigma);
325  if (gamma <= MagickEpsilon)
326  return(3UL);
327  alpha=PerceptibleReciprocal(2.0*gamma*gamma);
328  beta=(double) PerceptibleReciprocal((double) Magick2PI*gamma*gamma);
329  for (width=5; ; )
330  {
331  normalize=0.0;
332  j=(ssize_t) (width-1)/2;
333  for (v=(-j); v <= j; v++)
334  for (u=(-j); u <= j; u++)
335  normalize+=exp(-((double) (u*u+v*v))*alpha)*beta;
336  value=exp(-((double) (j*j))*alpha)*beta/normalize;
337  if ((value < QuantumScale) || (value < MagickEpsilon))
338  break;
339  width+=2;
340  }
341  return((size_t) (width-2));
342 }
343 
344 MagickPrivate size_t GetOptimalKernelWidth(const double radius,
345  const double sigma)
346 {
347  return(GetOptimalKernelWidth1D(radius,sigma));
348 }
_AffineMatrix
Definition: geometry.h:94
_RandomInfo
Definition: random.c:83