/src/imagemagick/MagickCore/gem.c

Source
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%                              GGGG  EEEEE  M   M                             %
%                             G      E      MM MM                             %
%                             G GG   EEE    M M M                             %
%                             G   G  E      M   M                             %
%                              GGGG  EEEEE  M   M                             %
%                                                                             %
%                                                                             %
%                    Graphic Gems - Graphic Support Methods                   %
%                                                                             %
%                               Software Design                               %
%                                    Cristy                                   %
%                                 August 1996                                 %
%                                                                             %
%                                                                             %
%  Copyright @ 1999 ImageMagick Studio LLC, a non-profit organization         %
%  dedicated to making software imaging solutions freely available.           %
%                                                                             %
%  You may not use this file except in compliance with the License.  You may  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    https://imagemagick.org/license/                                         %
%                                                                             %
%  Unless required by applicable law or agreed to in writing, software        %
%  distributed under the License is distributed on an "AS IS" BASIS,          %
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
%  See the License for the specific language governing permissions and        %
%  limitations under the License.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
*/

/*
  Include declarations.
*/
#include "MagickCore/studio.h"
#include "MagickCore/color-private.h"
#include "MagickCore/draw.h"
#include "MagickCore/gem.h"
#include "MagickCore/gem-private.h"
#include "MagickCore/image.h"
#include "MagickCore/image-private.h"
#include "MagickCore/log.h"
#include "MagickCore/memory_.h"
#include "MagickCore/pixel-accessor.h"
#include "MagickCore/quantum.h"
#include "MagickCore/quantum-private.h"
#include "MagickCore/random_.h"
#include "MagickCore/resize.h"
#include "MagickCore/transform.h"
#include "MagickCore/signature-private.h"

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   E x p a n d A f f i n e                                                   %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ExpandAffine() computes the affine's expansion factor, i.e. the square root
%  of the factor by which the affine transform affects area. In an affine
%  transform composed of scaling, rotation, shearing, and translation, returns
%  the amount of scaling.
%
%  The format of the ExpandAffine method is:
%
%      double ExpandAffine(const AffineMatrix *affine)
%
%  A description of each parameter follows:
%
%    o expansion: ExpandAffine returns the affine's expansion factor.
%
%    o affine: A pointer the affine transform of type AffineMatrix.
%
*/
MagickExport double ExpandAffine(const AffineMatrix *affine)
{
  assert(affine != (const AffineMatrix *) NULL);
  return(sqrt(fabs(affine->sx*affine->sy-affine->rx*affine->ry)));
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   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                         %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GenerateDifferentialNoise() generates differential noise.
%
%  The format of the GenerateDifferentialNoise method is:
%
%      double GenerateDifferentialNoise(RandomInfo *random_info,
%        const Quantum pixel,const NoiseType noise_type,const double attenuate)
%
%  A description of each parameter follows:
%
%    o random_info: the random info.
%
%    o pixel: noise is relative to this pixel value.
%
%    o noise_type: the type of noise.
%
%    o attenuate:  attenuate the noise.
%
*/
MagickPrivate double GenerateDifferentialNoise(RandomInfo *random_info,
  const Quantum pixel,const NoiseType noise_type,const double attenuate)
{
#define SigmaUniform  (attenuate*0.015625)
#define SigmaGaussian  (attenuate*0.015625)
#define SigmaImpulse  (attenuate*0.1)
#define SigmaLaplacian (attenuate*0.0390625)
#define SigmaMultiplicativeGaussian  (attenuate*0.5)
#define SigmaPoisson  (attenuate*12.5)
#define SigmaRandom  (attenuate)
#define TauGaussian  (attenuate*0.078125)

  double
    alpha,
    beta,
    noise,
    sigma;

  alpha=GetPseudoRandomValue(random_info);
  switch (noise_type)
  {
    case UniformNoise:
    default:
    {
      noise=(double) pixel+(double) QuantumRange*SigmaUniform*(alpha-0.5);
      break;
    }
    case GaussianNoise:
    {
      double
        gamma,
        tau;

      if (fabs(alpha) < MagickEpsilon)
        alpha=1.0;
      beta=GetPseudoRandomValue(random_info);
      gamma=sqrt(-2.0*log(alpha));
      sigma=gamma*cos((double) (2.0*MagickPI*beta));
      tau=gamma*sin((double) (2.0*MagickPI*beta));
      noise=(double) pixel+sqrt((double) pixel)*SigmaGaussian*sigma+
        (double) QuantumRange*TauGaussian*tau;
      break;
    }
    case ImpulseNoise:
    {
      if (alpha < (SigmaImpulse/2.0))
        noise=0.0;
      else
        if (alpha >= (1.0-(SigmaImpulse/2.0)))
          noise=(double) QuantumRange;
        else
          noise=(double) pixel;
      break;
    }
    case LaplacianNoise:
    {
      if (alpha <= 0.5)
        {
          if (alpha <= MagickEpsilon)
            noise=(double) (pixel-QuantumRange);
          else
            noise=(double) pixel+(double) QuantumRange*SigmaLaplacian*
              log(2.0*alpha)+0.5;
          break;
        }
      beta=1.0-alpha;
      if (beta <= (0.5*MagickEpsilon))
        noise=(double) (pixel+QuantumRange);
      else
        noise=(double) pixel-(double) QuantumRange*SigmaLaplacian*
          log(2.0*beta)+0.5;
      break;
    }
    case MultiplicativeGaussianNoise:
    {
      sigma=1.0;
      if (alpha > MagickEpsilon)
        sigma=sqrt(-2.0*log(alpha));
      beta=GetPseudoRandomValue(random_info);
      noise=(double) pixel+(double) pixel*SigmaMultiplicativeGaussian*sigma*
        cos((double) (2.0*MagickPI*beta))/2.0;
      break;
    }
    case PoissonNoise:
    {
      double
        poisson;

      ssize_t
        i;

      poisson=exp(-SigmaPoisson*QuantumScale*(double) pixel);
      for (i=0; alpha > poisson; i++)
      {
        beta=GetPseudoRandomValue(random_info);
        alpha*=beta;
      }
      noise=(double) QuantumRange*i*MagickSafeReciprocal(SigmaPoisson);
      break;
    }
    case RandomNoise:
    {
      noise=(double) QuantumRange*SigmaRandom*alpha;
      break;
    }
  }
  return(noise);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   G e t O p t i m a l K e r n e l W i d t h                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetOptimalKernelWidth() computes the optimal kernel radius for a convolution
%  filter.  Start with the minimum value of 3 pixels and walk out until we drop
%  below the threshold of one pixel numerical accuracy.
%
%  The format of the GetOptimalKernelWidth method is:
%
%      size_t GetOptimalKernelWidth(const double radius,
%        const double sigma)
%
%  A description of each parameter follows:
%
%    o width: GetOptimalKernelWidth returns the optimal width of a
%      convolution kernel.
%
%    o radius: the radius of the Gaussian, in pixels, not counting the center
%      pixel.
%
%    o sigma: the standard deviation of the Gaussian, in pixels.
%
*/
MagickPrivate size_t GetOptimalKernelWidth1D(const double radius,
  const double sigma)
{
  double
    alpha,
    beta,
    gamma,
    normalize,
    value;

  size_t
    width;

  ssize_t
    i,
    j;
 
  if (IsEventLogging() != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
  if (radius > MagickEpsilon)
    return((size_t) (2.0*ceil(radius)+1.0));
  gamma=fabs(sigma);
  if (gamma <= MagickEpsilon)
    return(3UL);
  alpha=MagickSafeReciprocal(2.0*gamma*gamma);
  beta=(double) MagickSafeReciprocal((double) MagickSQ2PI*gamma);
  for (width=5; ; )
  {
    normalize=0.0;
    j=(ssize_t) (width-1)/2;
    for (i=(-j); i <= j; i++)
      normalize+=exp(-((double) (i*i))*alpha)*beta;
    value=exp(-((double) (j*j))*alpha)*beta/normalize;
    if ((value < QuantumScale) || (value < MagickEpsilon))
      break;
    width+=2;
  }
  return((size_t) (width-2));
}

MagickPrivate size_t GetOptimalKernelWidth2D(const double radius,
  const double sigma)
{
  double
    alpha,
    beta,
    gamma,
    normalize,
    value;

  size_t
    width;

  ssize_t
    j,
    u,
    v;

  if (IsEventLogging() != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
  if (radius > MagickEpsilon)
    return((size_t) (2.0*ceil(radius)+1.0));
  gamma=fabs(sigma);
  if (gamma <= MagickEpsilon)
    return(3UL);
  alpha=MagickSafeReciprocal(2.0*gamma*gamma);
  beta=(double) MagickSafeReciprocal((double) Magick2PI*gamma*gamma);
  for (width=5; ; )
  {
    normalize=0.0;
    j=(ssize_t) (width-1)/2;
    for (v=(-j); v <= j; v++)
      for (u=(-j); u <= j; u++)
        normalize+=exp(-((double) (u*u+v*v))*alpha)*beta;
    value=exp(-((double) (j*j))*alpha)*beta/normalize;
    if ((value < QuantumScale) || (value < MagickEpsilon))
      break;
    width+=2;
  }
  return((size_t) (width-2));
}

MagickPrivate size_t  GetOptimalKernelWidth(const double radius,
  const double sigma)
{
  return(GetOptimalKernelWidth1D(radius,sigma));
}

Coverage Report

Created: 2026-03-31 06:56

Line	Count	Source
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/license/ %
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	2.04M	{
89	2.04M	assert(affine != (const AffineMatrix *) NULL);
90	2.04M	return(sqrt(fabs(affine->sxaffine->sy-affine->rxaffine->ry)));
91	2.04M	}
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	0	{
125	0	#define SigmaUniform (attenuate*0.015625)
126	0	#define SigmaGaussian (attenuate*0.015625)
127	0	#define SigmaImpulse (attenuate*0.1)
128	0	#define SigmaLaplacian (attenuate*0.0390625)
129	0	#define SigmaMultiplicativeGaussian (attenuate*0.5)
130	0	#define SigmaPoisson (attenuate*12.5)
131	0	#define SigmaRandom (attenuate)
132	0	#define TauGaussian (attenuate*0.078125)
133
134	0	double
135	0	alpha,
136	0	beta,
137	0	noise,
138	0	sigma;
139
140	0	alpha=GetPseudoRandomValue(random_info);
141	0	switch (noise_type)
142	0	{
143	0	case UniformNoise:
144	0	default:
145	0	{
146	0	noise=(double) pixel+(double) QuantumRangeSigmaUniform(alpha-0.5);
147	0	break;
148	0	}
149	0	case GaussianNoise:
150	0	{
151	0	double
152	0	gamma,
153	0	tau;
154
155	0	if (fabs(alpha) < MagickEpsilon)
156	0	alpha=1.0;
157	0	beta=GetPseudoRandomValue(random_info);
158	0	gamma=sqrt(-2.0*log(alpha));
159	0	sigma=gammacos((double) (2.0MagickPI*beta));
160	0	tau=gammasin((double) (2.0MagickPI*beta));
161	0	noise=(double) pixel+sqrt((double) pixel)SigmaGaussiansigma+
162	0	(double) QuantumRangeTauGaussiantau;
163	0	break;
164	0	}
165	0	case ImpulseNoise:
166	0	{
167	0	if (alpha < (SigmaImpulse/2.0))
168	0	noise=0.0;
169	0	else
170	0	if (alpha >= (1.0-(SigmaImpulse/2.0)))
171	0	noise=(double) QuantumRange;
172	0	else
173	0	noise=(double) pixel;
174	0	break;
175	0	}
176	0	case LaplacianNoise:
177	0	{
178	0	if (alpha <= 0.5)
179	0	{
180	0	if (alpha <= MagickEpsilon)
181	0	noise=(double) (pixel-QuantumRange);
182	0	else
183	0	noise=(double) pixel+(double) QuantumRangeSigmaLaplacian
184	0	log(2.0*alpha)+0.5;
185	0	break;
186	0	}
187	0	beta=1.0-alpha;
188	0	if (beta <= (0.5*MagickEpsilon))
189	0	noise=(double) (pixel+QuantumRange);
190	0	else
191	0	noise=(double) pixel-(double) QuantumRangeSigmaLaplacian
192	0	log(2.0*beta)+0.5;
193	0	break;
194	0	}
195	0	case MultiplicativeGaussianNoise:
196	0	{
197	0	sigma=1.0;
198	0	if (alpha > MagickEpsilon)
199	0	sigma=sqrt(-2.0*log(alpha));
200	0	beta=GetPseudoRandomValue(random_info);
201	0	noise=(double) pixel+(double) pixelSigmaMultiplicativeGaussiansigma*
202	0	cos((double) (2.0MagickPIbeta))/2.0;
203	0	break;
204	0	}
205	0	case PoissonNoise:
206	0	{
207	0	double
208	0	poisson;
209
210	0	ssize_t
211	0	i;
212
213	0	poisson=exp(-SigmaPoissonQuantumScale(double) pixel);
214	0	for (i=0; alpha > poisson; i++)
215	0	{
216	0	beta=GetPseudoRandomValue(random_info);
217	0	alpha*=beta;
218	0	}
219	0	noise=(double) QuantumRangeiMagickSafeReciprocal(SigmaPoisson);
220	0	break;
221	0	}
222	0	case RandomNoise:
223	0	{
224	0	noise=(double) QuantumRangeSigmaRandomalpha;
225	0	break;
226	0	}
227	0	}
228	0	return(noise);
229	0	}
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	0	{
265	0	double
266	0	alpha,
267	0	beta,
268	0	gamma,
269	0	normalize,
270	0	value;
271
272	0	size_t
273	0	width;
274
275	0	ssize_t
276	0	i,
277	0	j;
278
279	0	if (IsEventLogging() != MagickFalse)
280	0	(void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
281	0	if (radius > MagickEpsilon)
282	0	return((size_t) (2.0*ceil(radius)+1.0));
283	0	gamma=fabs(sigma);
284	0	if (gamma <= MagickEpsilon)
285	0	return(3UL);
286	0	alpha=MagickSafeReciprocal(2.0gammagamma);
287	0	beta=(double) MagickSafeReciprocal((double) MagickSQ2PI*gamma);
288	0	for (width=5; ; )
289	0	{
290	0	normalize=0.0;
291	0	j=(ssize_t) (width-1)/2;
292	0	for (i=(-j); i <= j; i++)
293	0	normalize+=exp(-((double) (ii))alpha)*beta;
294	0	value=exp(-((double) (jj))alpha)*beta/normalize;
295	0	if ((value < QuantumScale) \|\| (value < MagickEpsilon))
296	0	break;
297	0	width+=2;
298	0	}
299	0	return((size_t) (width-2));
300	0	}
301
302		MagickPrivate size_t GetOptimalKernelWidth2D(const double radius,
303		const double sigma)
304	0	{
305	0	double
306	0	alpha,
307	0	beta,
308	0	gamma,
309	0	normalize,
310	0	value;
311
312	0	size_t
313	0	width;
314
315	0	ssize_t
316	0	j,
317	0	u,
318	0	v;
319
320	0	if (IsEventLogging() != MagickFalse)
321	0	(void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
322	0	if (radius > MagickEpsilon)
323	0	return((size_t) (2.0*ceil(radius)+1.0));
324	0	gamma=fabs(sigma);
325	0	if (gamma <= MagickEpsilon)
326	0	return(3UL);
327	0	alpha=MagickSafeReciprocal(2.0gammagamma);
328	0	beta=(double) MagickSafeReciprocal((double) Magick2PIgammagamma);
329	0	for (width=5; ; )
330	0	{
331	0	normalize=0.0;
332	0	j=(ssize_t) (width-1)/2;
333	0	for (v=(-j); v <= j; v++)
334	0	for (u=(-j); u <= j; u++)
335	0	normalize+=exp(-((double) (uu+vv))alpha)beta;
336	0	value=exp(-((double) (jj))alpha)*beta/normalize;
337	0	if ((value < QuantumScale) \|\| (value < MagickEpsilon))
338	0	break;
339	0	width+=2;
340	0	}
341	0	return((size_t) (width-2));
342	0	}
343
344		MagickPrivate size_t GetOptimalKernelWidth(const double radius,
345		const double sigma)
346	0	{
347	0	return(GetOptimalKernelWidth1D(radius,sigma));
348	0	}