MagickCore  7.1.1-43
Convert, Edit, Or Compose Bitmap Images
threshold.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % TTTTT H H RRRR EEEEE SSSSS H H OOO L DDDD %
7 % T H H R R E SS H H O O L D D %
8 % T HHHHH RRRR EEE SSS HHHHH O O L D D %
9 % T H H R R E SS H H O O L D D %
10 % T H H R R EEEEE SSSSS H H OOO LLLLL DDDD %
11 % %
12 % %
13 % MagickCore Image Threshold Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % October 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/artifact.h"
45 #include "MagickCore/blob.h"
46 #include "MagickCore/cache-view.h"
47 #include "MagickCore/color.h"
48 #include "MagickCore/color-private.h"
49 #include "MagickCore/colormap.h"
50 #include "MagickCore/colorspace.h"
51 #include "MagickCore/colorspace-private.h"
52 #include "MagickCore/configure.h"
53 #include "MagickCore/constitute.h"
54 #include "MagickCore/decorate.h"
55 #include "MagickCore/draw.h"
56 #include "MagickCore/enhance.h"
57 #include "MagickCore/exception.h"
58 #include "MagickCore/exception-private.h"
59 #include "MagickCore/effect.h"
60 #include "MagickCore/fx.h"
61 #include "MagickCore/gem.h"
62 #include "MagickCore/gem-private.h"
63 #include "MagickCore/geometry.h"
64 #include "MagickCore/image-private.h"
65 #include "MagickCore/list.h"
66 #include "MagickCore/log.h"
67 #include "MagickCore/memory_.h"
68 #include "MagickCore/monitor.h"
69 #include "MagickCore/monitor-private.h"
70 #include "MagickCore/montage.h"
71 #include "MagickCore/option.h"
72 #include "MagickCore/pixel-accessor.h"
73 #include "MagickCore/property.h"
74 #include "MagickCore/quantize.h"
75 #include "MagickCore/quantum.h"
76 #include "MagickCore/quantum-private.h"
77 #include "MagickCore/random_.h"
78 #include "MagickCore/random-private.h"
79 #include "MagickCore/resize.h"
80 #include "MagickCore/resource_.h"
81 #include "MagickCore/segment.h"
82 #include "MagickCore/shear.h"
83 #include "MagickCore/signature-private.h"
84 #include "MagickCore/string_.h"
85 #include "MagickCore/string-private.h"
86 #include "MagickCore/thread-private.h"
87 #include "MagickCore/threshold.h"
88 #include "MagickCore/token.h"
89 #include "MagickCore/transform.h"
90 #include "MagickCore/xml-tree.h"
91 #include "MagickCore/xml-tree-private.h"
92 
93 /*
94  Define declarations.
95 */
96 #define ThresholdsFilename "thresholds.xml"
97 
98 /*
99  Typedef declarations.
100 */
102 {
103  char
104  *map_id,
105  *description;
106 
107  size_t
108  width,
109  height;
110 
111  ssize_t
112  divisor,
113  *levels;
114 };
115 
116 /*
117  Static declarations.
118 */
119 #if MAGICKCORE_ZERO_CONFIGURATION_SUPPORT
120  #include "MagickCore/threshold-map.h"
121 #else
122 static const char *const
123  BuiltinMap=
124  "<?xml version=\"1.0\"?>"
125  "<thresholds>"
126  " <threshold map=\"threshold\" alias=\"1x1\">"
127  " <description>Threshold 1x1 (non-dither)</description>"
128  " <levels width=\"1\" height=\"1\" divisor=\"2\">"
129  " 1"
130  " </levels>"
131  " </threshold>"
132  " <threshold map=\"checks\" alias=\"2x1\">"
133  " <description>Checkerboard 2x1 (dither)</description>"
134  " <levels width=\"2\" height=\"2\" divisor=\"3\">"
135  " 1 2"
136  " 2 1"
137  " </levels>"
138  " </threshold>"
139  "</thresholds>";
140 #endif
141 
142 /*
143  Forward declarations.
144 */
145 static ThresholdMap
146  *GetThresholdMapFile(const char *,const char *,const char *,ExceptionInfo *);
147 
148 /*
149 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
150 % %
151 % %
152 % %
153 % A d a p t i v e T h r e s h o l d I m a g e %
154 % %
155 % %
156 % %
157 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
158 %
159 % AdaptiveThresholdImage() selects an individual threshold for each pixel
160 % based on the range of intensity values in its local neighborhood. This
161 % allows for thresholding of an image whose global intensity histogram
162 % doesn't contain distinctive peaks.
163 %
164 % The format of the AdaptiveThresholdImage method is:
165 %
166 % Image *AdaptiveThresholdImage(const Image *image,const size_t width,
167 % const size_t height,const double bias,ExceptionInfo *exception)
168 %
169 % A description of each parameter follows:
170 %
171 % o image: the image.
172 %
173 % o width: the width of the local neighborhood.
174 %
175 % o height: the height of the local neighborhood.
176 %
177 % o bias: the mean bias.
178 %
179 % o exception: return any errors or warnings in this structure.
180 %
181 */
182 MagickExport Image *AdaptiveThresholdImage(const Image *image,
183  const size_t width,const size_t height,const double bias,
184  ExceptionInfo *exception)
185 {
186 #define AdaptiveThresholdImageTag "AdaptiveThreshold/Image"
187 
188  CacheView
189  *image_view,
190  *threshold_view;
191 
192  Image
193  *threshold_image;
194 
195  MagickBooleanType
196  status;
197 
198  MagickOffsetType
199  progress;
200 
201  MagickSizeType
202  number_pixels;
203 
204  ssize_t
205  y;
206 
207  /*
208  Initialize threshold image attributes.
209  */
210  assert(image != (Image *) NULL);
211  assert(image->signature == MagickCoreSignature);
212  assert(exception != (ExceptionInfo *) NULL);
213  assert(exception->signature == MagickCoreSignature);
214  if (IsEventLogging() != MagickFalse)
215  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
216  threshold_image=CloneImage(image,0,0,MagickTrue,exception);
217  if (threshold_image == (Image *) NULL)
218  return((Image *) NULL);
219  if ((width == 0) || (height == 0))
220  return(threshold_image);
221  status=SetImageStorageClass(threshold_image,DirectClass,exception);
222  if (status == MagickFalse)
223  {
224  threshold_image=DestroyImage(threshold_image);
225  return((Image *) NULL);
226  }
227  /*
228  Threshold image.
229  */
230  status=MagickTrue;
231  progress=0;
232  number_pixels=(MagickSizeType) width*height;
233  image_view=AcquireVirtualCacheView(image,exception);
234  threshold_view=AcquireAuthenticCacheView(threshold_image,exception);
235 #if defined(MAGICKCORE_OPENMP_SUPPORT)
236  #pragma omp parallel for schedule(static) shared(progress,status) \
237  magick_number_threads(image,threshold_image,image->rows,1)
238 #endif
239  for (y=0; y < (ssize_t) image->rows; y++)
240  {
241  double
242  channel_bias[MaxPixelChannels],
243  channel_sum[MaxPixelChannels];
244 
245  const Quantum
246  *magick_restrict p,
247  *magick_restrict pixels;
248 
249  Quantum
250  *magick_restrict q;
251 
252  ssize_t
253  center,
254  i,
255  u,
256  v,
257  x;
258 
259  if (status == MagickFalse)
260  continue;
261  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
262  (height/2L),image->columns+width,height,exception);
263  q=QueueCacheViewAuthenticPixels(threshold_view,0,y,threshold_image->columns,
264  1,exception);
265  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
266  {
267  status=MagickFalse;
268  continue;
269  }
270  center=(ssize_t) GetPixelChannels(image)*((ssize_t) image->columns+
271  (ssize_t) width)*((ssize_t) height/2L)+(ssize_t) GetPixelChannels(image)*
272  ((ssize_t) width/2);
273  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
274  {
275  PixelChannel channel = GetPixelChannelChannel(image,i);
276  PixelTrait traits = GetPixelChannelTraits(image,channel);
277  PixelTrait threshold_traits=GetPixelChannelTraits(threshold_image,
278  channel);
279  if ((traits == UndefinedPixelTrait) ||
280  (threshold_traits == UndefinedPixelTrait))
281  continue;
282  if ((threshold_traits & CopyPixelTrait) != 0)
283  {
284  SetPixelChannel(threshold_image,channel,p[center+i],q);
285  continue;
286  }
287  pixels=p;
288  channel_bias[channel]=0.0;
289  channel_sum[channel]=0.0;
290  for (v=0; v < (ssize_t) height; v++)
291  {
292  for (u=0; u < (ssize_t) width; u++)
293  {
294  if (u == (ssize_t) (width-1))
295  channel_bias[channel]+=(double) pixels[i];
296  channel_sum[channel]+=(double) pixels[i];
297  pixels+=(ptrdiff_t) GetPixelChannels(image);
298  }
299  pixels+=(ptrdiff_t) GetPixelChannels(image)*image->columns;
300  }
301  }
302  for (x=0; x < (ssize_t) image->columns; x++)
303  {
304  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
305  {
306  double
307  mean;
308 
309  PixelChannel channel = GetPixelChannelChannel(image,i);
310  PixelTrait traits = GetPixelChannelTraits(image,channel);
311  PixelTrait threshold_traits=GetPixelChannelTraits(threshold_image,
312  channel);
313  if ((traits == UndefinedPixelTrait) ||
314  (threshold_traits == UndefinedPixelTrait))
315  continue;
316  if ((threshold_traits & CopyPixelTrait) != 0)
317  {
318  SetPixelChannel(threshold_image,channel,p[center+i],q);
319  continue;
320  }
321  channel_sum[channel]-=channel_bias[channel];
322  channel_bias[channel]=0.0;
323  pixels=p;
324  for (v=0; v < (ssize_t) height; v++)
325  {
326  channel_bias[channel]+=(double) pixels[i];
327  pixels+=(width-1)*GetPixelChannels(image);
328  channel_sum[channel]+=(double) pixels[i];
329  pixels+=(ptrdiff_t) GetPixelChannels(image)*(image->columns+1);
330  }
331  mean=(double) (channel_sum[channel]/number_pixels+bias);
332  SetPixelChannel(threshold_image,channel,(Quantum) ((double)
333  p[center+i] <= mean ? 0 : QuantumRange),q);
334  }
335  p+=(ptrdiff_t) GetPixelChannels(image);
336  q+=(ptrdiff_t) GetPixelChannels(threshold_image);
337  }
338  if (SyncCacheViewAuthenticPixels(threshold_view,exception) == MagickFalse)
339  status=MagickFalse;
340  if (image->progress_monitor != (MagickProgressMonitor) NULL)
341  {
342  MagickBooleanType
343  proceed;
344 
345 #if defined(MAGICKCORE_OPENMP_SUPPORT)
346  #pragma omp atomic
347 #endif
348  progress++;
349  proceed=SetImageProgress(image,AdaptiveThresholdImageTag,progress,
350  image->rows);
351  if (proceed == MagickFalse)
352  status=MagickFalse;
353  }
354  }
355  threshold_image->type=image->type;
356  threshold_view=DestroyCacheView(threshold_view);
357  image_view=DestroyCacheView(image_view);
358  if (status == MagickFalse)
359  threshold_image=DestroyImage(threshold_image);
360  return(threshold_image);
361 }
362 
363 /*
364 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
365 % %
366 % %
367 % %
368 % A u t o T h r e s h o l d I m a g e %
369 % %
370 % %
371 % %
372 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
373 %
374 % AutoThresholdImage() automatically performs image thresholding
375 % dependent on which method you specify.
376 %
377 % The format of the AutoThresholdImage method is:
378 %
379 % MagickBooleanType AutoThresholdImage(Image *image,
380 % const AutoThresholdMethod method,ExceptionInfo *exception)
381 %
382 % A description of each parameter follows:
383 %
384 % o image: The image to auto-threshold.
385 %
386 % o method: choose from Kapur, OTSU, or Triangle.
387 %
388 % o exception: return any errors or warnings in this structure.
389 %
390 */
391 
392 static double KapurThreshold(const Image *image,const double *histogram,
393  ExceptionInfo *exception)
394 {
395 #define MaxIntensity 255
396 
397  double
398  *black_entropy,
399  *cumulative_histogram,
400  entropy,
401  epsilon,
402  maximum_entropy,
403  *white_entropy;
404 
405  ssize_t
406  i,
407  j;
408 
409  size_t
410  threshold;
411 
412  /*
413  Compute optimal threshold from the entropy of the histogram.
414  */
415  cumulative_histogram=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
416  sizeof(*cumulative_histogram));
417  black_entropy=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
418  sizeof(*black_entropy));
419  white_entropy=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
420  sizeof(*white_entropy));
421  if ((cumulative_histogram == (double *) NULL) ||
422  (black_entropy == (double *) NULL) || (white_entropy == (double *) NULL))
423  {
424  if (white_entropy != (double *) NULL)
425  white_entropy=(double *) RelinquishMagickMemory(white_entropy);
426  if (black_entropy != (double *) NULL)
427  black_entropy=(double *) RelinquishMagickMemory(black_entropy);
428  if (cumulative_histogram != (double *) NULL)
429  cumulative_histogram=(double *)
430  RelinquishMagickMemory(cumulative_histogram);
431  (void) ThrowMagickException(exception,GetMagickModule(),
432  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
433  return(-1.0);
434  }
435  /*
436  Entropy for black and white parts of the histogram.
437  */
438  cumulative_histogram[0]=histogram[0];
439  for (i=1; i <= MaxIntensity; i++)
440  cumulative_histogram[i]=cumulative_histogram[i-1]+histogram[i];
441  epsilon=MagickMinimumValue;
442  for (j=0; j <= MaxIntensity; j++)
443  {
444  /*
445  Black entropy.
446  */
447  black_entropy[j]=0.0;
448  if (cumulative_histogram[j] > epsilon)
449  {
450  entropy=0.0;
451  for (i=0; i <= j; i++)
452  if (histogram[i] > epsilon)
453  entropy-=histogram[i]/cumulative_histogram[j]*
454  log(histogram[i]/cumulative_histogram[j]);
455  black_entropy[j]=entropy;
456  }
457  /*
458  White entropy.
459  */
460  white_entropy[j]=0.0;
461  if ((1.0-cumulative_histogram[j]) > epsilon)
462  {
463  entropy=0.0;
464  for (i=j+1; i <= MaxIntensity; i++)
465  if (histogram[i] > epsilon)
466  entropy-=histogram[i]/(1.0-cumulative_histogram[j])*
467  log(histogram[i]/(1.0-cumulative_histogram[j]));
468  white_entropy[j]=entropy;
469  }
470  }
471  /*
472  Find histogram bin with maximum entropy.
473  */
474  maximum_entropy=black_entropy[0]+white_entropy[0];
475  threshold=0;
476  for (j=1; j <= MaxIntensity; j++)
477  if ((black_entropy[j]+white_entropy[j]) > maximum_entropy)
478  {
479  maximum_entropy=black_entropy[j]+white_entropy[j];
480  threshold=(size_t) j;
481  }
482  /*
483  Free resources.
484  */
485  white_entropy=(double *) RelinquishMagickMemory(white_entropy);
486  black_entropy=(double *) RelinquishMagickMemory(black_entropy);
487  cumulative_histogram=(double *) RelinquishMagickMemory(cumulative_histogram);
488  return(100.0*threshold/MaxIntensity);
489 }
490 
491 static double OTSUThreshold(const Image *image,const double *histogram,
492  ExceptionInfo *exception)
493 {
494  double
495  max_sigma,
496  *myu,
497  *omega,
498  *probability,
499  *sigma,
500  threshold;
501 
502  ssize_t
503  i;
504 
505  /*
506  Compute optimal threshold from maximization of inter-class variance.
507  */
508  myu=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*myu));
509  omega=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*omega));
510  probability=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
511  sizeof(*probability));
512  sigma=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*sigma));
513  if ((myu == (double *) NULL) || (omega == (double *) NULL) ||
514  (probability == (double *) NULL) || (sigma == (double *) NULL))
515  {
516  if (sigma != (double *) NULL)
517  sigma=(double *) RelinquishMagickMemory(sigma);
518  if (probability != (double *) NULL)
519  probability=(double *) RelinquishMagickMemory(probability);
520  if (omega != (double *) NULL)
521  omega=(double *) RelinquishMagickMemory(omega);
522  if (myu != (double *) NULL)
523  myu=(double *) RelinquishMagickMemory(myu);
524  (void) ThrowMagickException(exception,GetMagickModule(),
525  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
526  return(-1.0);
527  }
528  /*
529  Calculate probability density.
530  */
531  for (i=0; i <= (ssize_t) MaxIntensity; i++)
532  probability[i]=histogram[i];
533  /*
534  Generate probability of graylevels and mean value for separation.
535  */
536  omega[0]=probability[0];
537  myu[0]=0.0;
538  for (i=1; i <= (ssize_t) MaxIntensity; i++)
539  {
540  omega[i]=omega[i-1]+probability[i];
541  myu[i]=myu[i-1]+i*probability[i];
542  }
543  /*
544  Sigma maximization: inter-class variance and compute optimal threshold.
545  */
546  threshold=0;
547  max_sigma=0.0;
548  for (i=0; i < (ssize_t) MaxIntensity; i++)
549  {
550  sigma[i]=0.0;
551  if ((omega[i] != 0.0) && (omega[i] != 1.0))
552  sigma[i]=pow(myu[MaxIntensity]*omega[i]-myu[i],2.0)/(omega[i]*(1.0-
553  omega[i]));
554  if (sigma[i] > max_sigma)
555  {
556  max_sigma=sigma[i];
557  threshold=(double) i;
558  }
559  }
560  /*
561  Free resources.
562  */
563  myu=(double *) RelinquishMagickMemory(myu);
564  omega=(double *) RelinquishMagickMemory(omega);
565  probability=(double *) RelinquishMagickMemory(probability);
566  sigma=(double *) RelinquishMagickMemory(sigma);
567  return(100.0*threshold/MaxIntensity);
568 }
569 
570 static double TriangleThreshold(const double *histogram)
571 {
572  double
573  a,
574  b,
575  c,
576  count,
577  distance,
578  inverse_ratio,
579  max_distance,
580  segment,
581  x1,
582  x2,
583  y1,
584  y2;
585 
586  ssize_t
587  i;
588 
589  ssize_t
590  end,
591  max,
592  start,
593  threshold;
594 
595  /*
596  Compute optimal threshold with triangle algorithm.
597  */
598  start=0; /* find start bin, first bin not zero count */
599  for (i=0; i <= (ssize_t) MaxIntensity; i++)
600  if (histogram[i] > 0.0)
601  {
602  start=i;
603  break;
604  }
605  end=0; /* find end bin, last bin not zero count */
606  for (i=(ssize_t) MaxIntensity; i >= 0; i--)
607  if (histogram[i] > 0.0)
608  {
609  end=i;
610  break;
611  }
612  max=0; /* find max bin, bin with largest count */
613  count=0.0;
614  for (i=0; i <= (ssize_t) MaxIntensity; i++)
615  if (histogram[i] > count)
616  {
617  max=i;
618  count=histogram[i];
619  }
620  /*
621  Compute threshold at split point.
622  */
623  x1=(double) max;
624  y1=histogram[max];
625  x2=(double) end;
626  if ((max-start) >= (end-max))
627  x2=(double) start;
628  y2=0.0;
629  a=y1-y2;
630  b=x2-x1;
631  c=(-1.0)*(a*x1+b*y1);
632  inverse_ratio=1.0/sqrt(a*a+b*b+c*c);
633  threshold=0;
634  max_distance=0.0;
635  if (x2 == (double) start)
636  for (i=start; i < max; i++)
637  {
638  segment=inverse_ratio*(a*i+b*histogram[i]+c);
639  distance=sqrt(segment*segment);
640  if ((distance > max_distance) && (segment > 0.0))
641  {
642  threshold=i;
643  max_distance=distance;
644  }
645  }
646  else
647  for (i=end; i > max; i--)
648  {
649  segment=inverse_ratio*(a*i+b*histogram[i]+c);
650  distance=sqrt(segment*segment);
651  if ((distance > max_distance) && (segment < 0.0))
652  {
653  threshold=i;
654  max_distance=distance;
655  }
656  }
657  return(100.0*threshold/MaxIntensity);
658 }
659 
660 MagickExport MagickBooleanType AutoThresholdImage(Image *image,
661  const AutoThresholdMethod method,ExceptionInfo *exception)
662 {
663  CacheView
664  *image_view;
665 
666  char
667  property[MagickPathExtent];
668 
669  double
670  gamma,
671  *histogram,
672  sum,
673  threshold;
674 
675  MagickBooleanType
676  status;
677 
678  ssize_t
679  i;
680 
681  ssize_t
682  y;
683 
684  /*
685  Form histogram.
686  */
687  assert(image != (Image *) NULL);
688  assert(image->signature == MagickCoreSignature);
689  if (IsEventLogging() != MagickFalse)
690  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
691  histogram=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
692  sizeof(*histogram));
693  if (histogram == (double *) NULL)
694  ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
695  image->filename);
696  status=MagickTrue;
697  (void) memset(histogram,0,(MaxIntensity+1UL)*sizeof(*histogram));
698  image_view=AcquireVirtualCacheView(image,exception);
699  for (y=0; y < (ssize_t) image->rows; y++)
700  {
701  const Quantum
702  *magick_restrict p;
703 
704  ssize_t
705  x;
706 
707  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
708  if (p == (const Quantum *) NULL)
709  break;
710  for (x=0; x < (ssize_t) image->columns; x++)
711  {
712  double intensity = GetPixelIntensity(image,p);
713  histogram[ScaleQuantumToChar(ClampToQuantum(intensity))]++;
714  p+=(ptrdiff_t) GetPixelChannels(image);
715  }
716  }
717  image_view=DestroyCacheView(image_view);
718  /*
719  Normalize histogram.
720  */
721  sum=0.0;
722  for (i=0; i <= (ssize_t) MaxIntensity; i++)
723  sum+=histogram[i];
724  gamma=PerceptibleReciprocal(sum);
725  for (i=0; i <= (ssize_t) MaxIntensity; i++)
726  histogram[i]=gamma*histogram[i];
727  /*
728  Discover threshold from histogram.
729  */
730  switch (method)
731  {
732  case KapurThresholdMethod:
733  {
734  threshold=KapurThreshold(image,histogram,exception);
735  break;
736  }
737  case OTSUThresholdMethod:
738  default:
739  {
740  threshold=OTSUThreshold(image,histogram,exception);
741  break;
742  }
743  case TriangleThresholdMethod:
744  {
745  threshold=TriangleThreshold(histogram);
746  break;
747  }
748  }
749  histogram=(double *) RelinquishMagickMemory(histogram);
750  if (threshold < 0.0)
751  status=MagickFalse;
752  if (status == MagickFalse)
753  return(MagickFalse);
754  /*
755  Threshold image.
756  */
757  (void) FormatLocaleString(property,MagickPathExtent,"%g%%",threshold);
758  (void) SetImageProperty(image,"auto-threshold:threshold",property,exception);
759  if (IsStringTrue(GetImageArtifact(image,"auto-threshold:verbose")) != MagickFalse)
760  (void) FormatLocaleFile(stdout,"%.*g%%\n",GetMagickPrecision(),threshold);
761  return(BilevelImage(image,(double) QuantumRange*threshold/100.0,exception));
762 }
763 
764 /*
765 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
766 % %
767 % %
768 % %
769 % B i l e v e l I m a g e %
770 % %
771 % %
772 % %
773 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
774 %
775 % BilevelImage() changes the value of individual pixels based on the
776 % intensity of each pixel channel. The result is a high-contrast image.
777 %
778 % More precisely each channel value of the image is 'thresholded' so that if
779 % it is equal to or less than the given value it is set to zero, while any
780 % value greater than that give is set to it maximum or QuantumRange.
781 %
782 % This function is what is used to implement the "-threshold" operator for
783 % the command line API.
784 %
785 % If the default channel setting is given the image is thresholded using just
786 % the gray 'intensity' of the image, rather than the individual channels.
787 %
788 % The format of the BilevelImage method is:
789 %
790 % MagickBooleanType BilevelImage(Image *image,const double threshold,
791 % ExceptionInfo *exception)
792 %
793 % A description of each parameter follows:
794 %
795 % o image: the image.
796 %
797 % o threshold: define the threshold values.
798 %
799 % o exception: return any errors or warnings in this structure.
800 %
801 % Aside: You can get the same results as operator using LevelImages()
802 % with the 'threshold' value for both the black_point and the white_point.
803 %
804 */
805 MagickExport MagickBooleanType BilevelImage(Image *image,const double threshold,
806  ExceptionInfo *exception)
807 {
808 #define ThresholdImageTag "Threshold/Image"
809 
810  CacheView
811  *image_view;
812 
813  MagickBooleanType
814  status;
815 
816  MagickOffsetType
817  progress;
818 
819  ssize_t
820  y;
821 
822  assert(image != (Image *) NULL);
823  assert(image->signature == MagickCoreSignature);
824  if (IsEventLogging() != MagickFalse)
825  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
826  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
827  return(MagickFalse);
828  if (IsGrayColorspace(image->colorspace) == MagickFalse)
829  (void) SetImageColorspace(image,sRGBColorspace,exception);
830  /*
831  Bilevel threshold image.
832  */
833  status=MagickTrue;
834  progress=0;
835  image_view=AcquireAuthenticCacheView(image,exception);
836 #if defined(MAGICKCORE_OPENMP_SUPPORT)
837  #pragma omp parallel for schedule(static) shared(progress,status) \
838  magick_number_threads(image,image,image->rows,2)
839 #endif
840  for (y=0; y < (ssize_t) image->rows; y++)
841  {
842  ssize_t
843  x;
844 
845  Quantum
846  *magick_restrict q;
847 
848  if (status == MagickFalse)
849  continue;
850  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
851  if (q == (Quantum *) NULL)
852  {
853  status=MagickFalse;
854  continue;
855  }
856  for (x=0; x < (ssize_t) image->columns; x++)
857  {
858  double
859  pixel;
860 
861  ssize_t
862  i;
863 
864  pixel=GetPixelIntensity(image,q);
865  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
866  {
867  PixelChannel channel = GetPixelChannelChannel(image,i);
868  PixelTrait traits = GetPixelChannelTraits(image,channel);
869  if ((traits & UpdatePixelTrait) == 0)
870  continue;
871  if (image->channel_mask != AllChannels)
872  pixel=(double) q[i];
873  q[i]=(Quantum) (pixel <= threshold ? 0 : QuantumRange);
874  }
875  q+=(ptrdiff_t) GetPixelChannels(image);
876  }
877  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
878  status=MagickFalse;
879  if (image->progress_monitor != (MagickProgressMonitor) NULL)
880  {
881  MagickBooleanType
882  proceed;
883 
884 #if defined(MAGICKCORE_OPENMP_SUPPORT)
885  #pragma omp atomic
886 #endif
887  progress++;
888  proceed=SetImageProgress(image,ThresholdImageTag,progress++,
889  image->rows);
890  if (proceed == MagickFalse)
891  status=MagickFalse;
892  }
893  }
894  image_view=DestroyCacheView(image_view);
895  return(status);
896 }
897 
898 /*
899 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
900 % %
901 % %
902 % %
903 % B l a c k T h r e s h o l d I m a g e %
904 % %
905 % %
906 % %
907 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
908 %
909 % BlackThresholdImage() is like ThresholdImage() but forces all pixels below
910 % the threshold into black while leaving all pixels at or above the threshold
911 % unchanged.
912 %
913 % The format of the BlackThresholdImage method is:
914 %
915 % MagickBooleanType BlackThresholdImage(Image *image,
916 % const char *threshold,ExceptionInfo *exception)
917 %
918 % A description of each parameter follows:
919 %
920 % o image: the image.
921 %
922 % o threshold: define the threshold value.
923 %
924 % o exception: return any errors or warnings in this structure.
925 %
926 */
927 MagickExport MagickBooleanType BlackThresholdImage(Image *image,
928  const char *thresholds,ExceptionInfo *exception)
929 {
930 #define ThresholdImageTag "Threshold/Image"
931 
932  CacheView
933  *image_view;
934 
936  geometry_info;
937 
938  MagickBooleanType
939  status;
940 
941  MagickOffsetType
942  progress;
943 
944  PixelInfo
945  threshold;
946 
947  MagickStatusType
948  flags;
949 
950  ssize_t
951  y;
952 
953  assert(image != (Image *) NULL);
954  assert(image->signature == MagickCoreSignature);
955  if (IsEventLogging() != MagickFalse)
956  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
957  if (thresholds == (const char *) NULL)
958  return(MagickTrue);
959  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
960  return(MagickFalse);
961  if (IsGrayColorspace(image->colorspace) != MagickFalse)
962  (void) SetImageColorspace(image,sRGBColorspace,exception);
963  GetPixelInfo(image,&threshold);
964  flags=ParseGeometry(thresholds,&geometry_info);
965  threshold.red=geometry_info.rho;
966  threshold.green=geometry_info.rho;
967  threshold.blue=geometry_info.rho;
968  threshold.black=geometry_info.rho;
969  threshold.alpha=100.0;
970  if ((flags & SigmaValue) != 0)
971  threshold.green=geometry_info.sigma;
972  if ((flags & XiValue) != 0)
973  threshold.blue=geometry_info.xi;
974  if ((flags & PsiValue) != 0)
975  threshold.alpha=geometry_info.psi;
976  if (threshold.colorspace == CMYKColorspace)
977  {
978  if ((flags & PsiValue) != 0)
979  threshold.black=geometry_info.psi;
980  if ((flags & ChiValue) != 0)
981  threshold.alpha=geometry_info.chi;
982  }
983  if ((flags & PercentValue) != 0)
984  {
985  threshold.red*=((MagickRealType) QuantumRange/100.0);
986  threshold.green*=((MagickRealType) QuantumRange/100.0);
987  threshold.blue*=((MagickRealType) QuantumRange/100.0);
988  threshold.black*=((MagickRealType) QuantumRange/100.0);
989  threshold.alpha*=((MagickRealType) QuantumRange/100.0);
990  }
991  /*
992  White threshold image.
993  */
994  status=MagickTrue;
995  progress=0;
996  image_view=AcquireAuthenticCacheView(image,exception);
997 #if defined(MAGICKCORE_OPENMP_SUPPORT)
998  #pragma omp parallel for schedule(static) shared(progress,status) \
999  magick_number_threads(image,image,image->rows,2)
1000 #endif
1001  for (y=0; y < (ssize_t) image->rows; y++)
1002  {
1003  ssize_t
1004  x;
1005 
1006  Quantum
1007  *magick_restrict q;
1008 
1009  if (status == MagickFalse)
1010  continue;
1011  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1012  if (q == (Quantum *) NULL)
1013  {
1014  status=MagickFalse;
1015  continue;
1016  }
1017  for (x=0; x < (ssize_t) image->columns; x++)
1018  {
1019  double
1020  pixel;
1021 
1022  ssize_t
1023  i;
1024 
1025  pixel=GetPixelIntensity(image,q);
1026  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1027  {
1028  PixelChannel channel = GetPixelChannelChannel(image,i);
1029  PixelTrait traits = GetPixelChannelTraits(image,channel);
1030  if ((traits & UpdatePixelTrait) == 0)
1031  continue;
1032  if (image->channel_mask != AllChannels)
1033  pixel=(double) q[i];
1034  if (pixel < GetPixelInfoChannel(&threshold,channel))
1035  q[i]=(Quantum) 0;
1036  }
1037  q+=(ptrdiff_t) GetPixelChannels(image);
1038  }
1039  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1040  status=MagickFalse;
1041  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1042  {
1043  MagickBooleanType
1044  proceed;
1045 
1046 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1047  #pragma omp atomic
1048 #endif
1049  progress++;
1050  proceed=SetImageProgress(image,ThresholdImageTag,progress,
1051  image->rows);
1052  if (proceed == MagickFalse)
1053  status=MagickFalse;
1054  }
1055  }
1056  image_view=DestroyCacheView(image_view);
1057  return(status);
1058 }
1059 
1060 /*
1061 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1062 % %
1063 % %
1064 % %
1065 % C l a m p I m a g e %
1066 % %
1067 % %
1068 % %
1069 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1070 %
1071 % ClampImage() set each pixel whose value is below zero to zero and any the
1072 % pixel whose value is above the quantum range to the quantum range (e.g.
1073 % 65535) otherwise the pixel value remains unchanged.
1074 %
1075 % The format of the ClampImage method is:
1076 %
1077 % MagickBooleanType ClampImage(Image *image,ExceptionInfo *exception)
1078 %
1079 % A description of each parameter follows:
1080 %
1081 % o image: the image.
1082 %
1083 % o exception: return any errors or warnings in this structure.
1084 %
1085 */
1086 
1087 MagickExport MagickBooleanType ClampImage(Image *image,ExceptionInfo *exception)
1088 {
1089 #define ClampImageTag "Clamp/Image"
1090 
1091  CacheView
1092  *image_view;
1093 
1094  MagickBooleanType
1095  status;
1096 
1097  MagickOffsetType
1098  progress;
1099 
1100  ssize_t
1101  y;
1102 
1103  assert(image != (Image *) NULL);
1104  assert(image->signature == MagickCoreSignature);
1105  if (IsEventLogging() != MagickFalse)
1106  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1107  if (image->storage_class == PseudoClass)
1108  {
1109  ssize_t
1110  i;
1111 
1112  PixelInfo
1113  *magick_restrict q;
1114 
1115  q=image->colormap;
1116  for (i=0; i < (ssize_t) image->colors; i++)
1117  {
1118  q->red=(double) ClampPixel(q->red);
1119  q->green=(double) ClampPixel(q->green);
1120  q->blue=(double) ClampPixel(q->blue);
1121  q->alpha=(double) ClampPixel(q->alpha);
1122  q++;
1123  }
1124  return(SyncImage(image,exception));
1125  }
1126  /*
1127  Clamp image.
1128  */
1129  status=MagickTrue;
1130  progress=0;
1131  image_view=AcquireAuthenticCacheView(image,exception);
1132 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1133  #pragma omp parallel for schedule(static) shared(progress,status) \
1134  magick_number_threads(image,image,image->rows,2)
1135 #endif
1136  for (y=0; y < (ssize_t) image->rows; y++)
1137  {
1138  ssize_t
1139  x;
1140 
1141  Quantum
1142  *magick_restrict q;
1143 
1144  if (status == MagickFalse)
1145  continue;
1146  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1147  if (q == (Quantum *) NULL)
1148  {
1149  status=MagickFalse;
1150  continue;
1151  }
1152  for (x=0; x < (ssize_t) image->columns; x++)
1153  {
1154  ssize_t
1155  i;
1156 
1157  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1158  {
1159  PixelChannel channel = GetPixelChannelChannel(image,i);
1160  PixelTrait traits = GetPixelChannelTraits(image,channel);
1161  if ((traits & UpdatePixelTrait) == 0)
1162  continue;
1163  q[i]=ClampPixel((MagickRealType) q[i]);
1164  }
1165  q+=(ptrdiff_t) GetPixelChannels(image);
1166  }
1167  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1168  status=MagickFalse;
1169  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1170  {
1171  MagickBooleanType
1172  proceed;
1173 
1174 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1175  #pragma omp atomic
1176 #endif
1177  progress++;
1178  proceed=SetImageProgress(image,ClampImageTag,progress,image->rows);
1179  if (proceed == MagickFalse)
1180  status=MagickFalse;
1181  }
1182  }
1183  image_view=DestroyCacheView(image_view);
1184  return(status);
1185 }
1186 
1187 /*
1188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1189 % %
1190 % %
1191 % %
1192 % C o l o r T h r e s h o l d I m a g e %
1193 % %
1194 % %
1195 % %
1196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1197 %
1198 % ColorThresholdImage() forces all pixels in the color range to white
1199 % otherwise black.
1200 %
1201 % The format of the ColorThresholdImage method is:
1202 %
1203 % MagickBooleanType ColorThresholdImage(Image *image,
1204 % const PixelInfo *start_color,const PixelInfo *stop_color,
1205 % ExceptionInfo *exception)
1206 %
1207 % A description of each parameter follows:
1208 %
1209 % o image: the image.
1210 %
1211 % o start_color, stop_color: define the start and stop color range. Any
1212 % pixel within the range returns white otherwise black.
1213 %
1214 % o exception: return any errors or warnings in this structure.
1215 %
1216 */
1217 MagickExport MagickBooleanType ColorThresholdImage(Image *image,
1218  const PixelInfo *start_color,const PixelInfo *stop_color,
1219  ExceptionInfo *exception)
1220 {
1221 #define ThresholdImageTag "Threshold/Image"
1222 
1223  CacheView
1224  *image_view;
1225 
1226  const char
1227  *artifact;
1228 
1229  IlluminantType
1230  illuminant = D65Illuminant;
1231 
1232  MagickBooleanType
1233  status;
1234 
1235  MagickOffsetType
1236  progress;
1237 
1238  PixelInfo
1239  start,
1240  stop;
1241 
1242  ssize_t
1243  y;
1244 
1245  /*
1246  Color threshold image.
1247  */
1248  assert(image != (Image *) NULL);
1249  assert(image->signature == MagickCoreSignature);
1250  if (IsEventLogging() != MagickFalse)
1251  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1252  status=AcquireImageColormap(image,2,exception);
1253  if (status == MagickFalse)
1254  return(status);
1255  artifact=GetImageArtifact(image,"color:illuminant");
1256  if (artifact != (const char *) NULL)
1257  {
1258  ssize_t
1259  illuminant_type;
1260 
1261  illuminant_type=ParseCommandOption(MagickIlluminantOptions,MagickFalse,
1262  artifact);
1263  if (illuminant_type < 0)
1264  illuminant=UndefinedIlluminant;
1265  else
1266  illuminant=(IlluminantType) illuminant_type;
1267  }
1268  start=(*start_color);
1269  stop=(*stop_color);
1270  switch (image->colorspace)
1271  {
1272  case HCLColorspace:
1273  {
1274  ConvertRGBToHCL(start_color->red,start_color->green,start_color->blue,
1275  &start.red,&start.green,&start.blue);
1276  ConvertRGBToHCL(stop_color->red,stop_color->green,stop_color->blue,
1277  &stop.red,&stop.green,&stop.blue);
1278  break;
1279  }
1280  case HSBColorspace:
1281  {
1282  ConvertRGBToHSB(start_color->red,start_color->green,start_color->blue,
1283  &start.red,&start.green,&start.blue);
1284  ConvertRGBToHSB(stop_color->red,stop_color->green,stop_color->blue,
1285  &stop.red,&stop.green,&stop.blue);
1286  break;
1287  }
1288  case HSLColorspace:
1289  {
1290  ConvertRGBToHSL(start_color->red,start_color->green,start_color->blue,
1291  &start.red,&start.green,&start.blue);
1292  ConvertRGBToHSL(stop_color->red,stop_color->green,stop_color->blue,
1293  &stop.red,&stop.green,&stop.blue);
1294  break;
1295  }
1296  case HSVColorspace:
1297  {
1298  ConvertRGBToHSV(start_color->red,start_color->green,start_color->blue,
1299  &start.red,&start.green,&start.blue);
1300  ConvertRGBToHSV(stop_color->red,stop_color->green,stop_color->blue,
1301  &stop.red,&stop.green,&stop.blue);
1302  break;
1303  }
1304  case HWBColorspace:
1305  {
1306  ConvertRGBToHWB(start_color->red,start_color->green,start_color->blue,
1307  &start.red,&start.green,&start.blue);
1308  ConvertRGBToHWB(stop_color->red,stop_color->green,stop_color->blue,
1309  &stop.red,&stop.green,&stop.blue);
1310  break;
1311  }
1312  case LabColorspace:
1313  {
1314  ConvertRGBToLab(start_color->red,start_color->green,start_color->blue,
1315  illuminant,&start.red,&start.green,&start.blue);
1316  ConvertRGBToLab(stop_color->red,stop_color->green,stop_color->blue,
1317  illuminant,&stop.red,&stop.green,&stop.blue);
1318  break;
1319  }
1320  default:
1321  {
1322  start.red*=QuantumScale;
1323  start.green*=QuantumScale;
1324  start.blue*=QuantumScale;
1325  stop.red*=QuantumScale;
1326  stop.green*=QuantumScale;
1327  stop.blue*=QuantumScale;
1328  break;
1329  }
1330  }
1331  start.red*=(double) QuantumRange;
1332  start.green*=(double) QuantumRange;
1333  start.blue*=(double) QuantumRange;
1334  stop.red*=(double) QuantumRange;
1335  stop.green*=(double) QuantumRange;
1336  stop.blue*=(double) QuantumRange;
1337  progress=0;
1338  image_view=AcquireAuthenticCacheView(image,exception);
1339 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1340  #pragma omp parallel for schedule(static) shared(progress,status) \
1341  magick_number_threads(image,image,image->rows,2)
1342 #endif
1343  for (y=0; y < (ssize_t) image->rows; y++)
1344  {
1345  ssize_t
1346  x;
1347 
1348  Quantum
1349  *magick_restrict q;
1350 
1351  if (status == MagickFalse)
1352  continue;
1353  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1354  if (q == (Quantum *) NULL)
1355  {
1356  status=MagickFalse;
1357  continue;
1358  }
1359  for (x=0; x < (ssize_t) image->columns; x++)
1360  {
1361  MagickBooleanType
1362  foreground = MagickTrue;
1363 
1364  ssize_t
1365  i;
1366 
1367  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1368  {
1369  PixelChannel channel = GetPixelChannelChannel(image,i);
1370  PixelTrait traits = GetPixelChannelTraits(image,channel);
1371  if ((traits & UpdatePixelTrait) == 0)
1372  continue;
1373  if (((double) q[i] < GetPixelInfoChannel(&start,channel)) ||
1374  ((double) q[i] > GetPixelInfoChannel(&stop,channel)))
1375  foreground=MagickFalse;
1376  }
1377  SetPixelIndex(image,(Quantum) (foreground != MagickFalse ? 1 : 0),q);
1378  q+=(ptrdiff_t) GetPixelChannels(image);
1379  }
1380  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1381  status=MagickFalse;
1382  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1383  {
1384  MagickBooleanType
1385  proceed;
1386 
1387 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1388  #pragma omp atomic
1389 #endif
1390  progress++;
1391  proceed=SetImageProgress(image,ThresholdImageTag,progress,
1392  image->rows);
1393  if (proceed == MagickFalse)
1394  status=MagickFalse;
1395  }
1396  }
1397  image_view=DestroyCacheView(image_view);
1398  image->colorspace=sRGBColorspace;
1399  return(SyncImage(image,exception));
1400 }
1401 
1402 /*
1403 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1404 % %
1405 % %
1406 % %
1407 % D e s t r o y T h r e s h o l d M a p %
1408 % %
1409 % %
1410 % %
1411 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1412 %
1413 % DestroyThresholdMap() de-allocate the given ThresholdMap
1414 %
1415 % The format of the ListThresholdMaps method is:
1416 %
1417 % ThresholdMap *DestroyThresholdMap(Threshold *map)
1418 %
1419 % A description of each parameter follows.
1420 %
1421 % o map: Pointer to the Threshold map to destroy
1422 %
1423 */
1424 MagickExport ThresholdMap *DestroyThresholdMap(ThresholdMap *map)
1425 {
1426  assert(map != (ThresholdMap *) NULL);
1427  if (map->map_id != (char *) NULL)
1428  map->map_id=DestroyString(map->map_id);
1429  if (map->description != (char *) NULL)
1430  map->description=DestroyString(map->description);
1431  if (map->levels != (ssize_t *) NULL)
1432  map->levels=(ssize_t *) RelinquishMagickMemory(map->levels);
1433  map=(ThresholdMap *) RelinquishMagickMemory(map);
1434  return(map);
1435 }
1436 
1437 /*
1438 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1439 % %
1440 % %
1441 % %
1442 % G e t T h r e s h o l d M a p %
1443 % %
1444 % %
1445 % %
1446 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1447 %
1448 % GetThresholdMap() loads and searches one or more threshold map files for the
1449 % map matching the given name or alias.
1450 %
1451 % The format of the GetThresholdMap method is:
1452 %
1453 % ThresholdMap *GetThresholdMap(const char *map_id,
1454 % ExceptionInfo *exception)
1455 %
1456 % A description of each parameter follows.
1457 %
1458 % o map_id: ID of the map to look for.
1459 %
1460 % o exception: return any errors or warnings in this structure.
1461 %
1462 */
1463 MagickExport ThresholdMap *GetThresholdMap(const char *map_id,
1464  ExceptionInfo *exception)
1465 {
1466  ThresholdMap
1467  *map;
1468 
1469  map=GetThresholdMapFile(BuiltinMap,"built-in",map_id,exception);
1470  if (map != (ThresholdMap *) NULL)
1471  return(map);
1472 #if !MAGICKCORE_ZERO_CONFIGURATION_SUPPORT
1473  {
1474  const StringInfo
1475  *option;
1476 
1478  *options;
1479 
1480  options=GetConfigureOptions(ThresholdsFilename,exception);
1481  option=(const StringInfo *) GetNextValueInLinkedList(options);
1482  while (option != (const StringInfo *) NULL)
1483  {
1484  map=GetThresholdMapFile((const char *) GetStringInfoDatum(option),
1485  GetStringInfoPath(option),map_id,exception);
1486  if (map != (ThresholdMap *) NULL)
1487  break;
1488  option=(const StringInfo *) GetNextValueInLinkedList(options);
1489  }
1490  options=DestroyConfigureOptions(options);
1491  }
1492 #endif
1493  return(map);
1494 }
1495 
1496 /*
1497 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1498 % %
1499 % %
1500 % %
1501 + G e t T h r e s h o l d M a p F i l e %
1502 % %
1503 % %
1504 % %
1505 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1506 %
1507 % GetThresholdMapFile() look for a given threshold map name or alias in the
1508 % given XML file data, and return the allocated the map when found.
1509 %
1510 % The format of the ListThresholdMaps method is:
1511 %
1512 % ThresholdMap *GetThresholdMap(const char *xml,const char *filename,
1513 % const char *map_id,ExceptionInfo *exception)
1514 %
1515 % A description of each parameter follows.
1516 %
1517 % o xml: The threshold map list in XML format.
1518 %
1519 % o filename: The threshold map XML filename.
1520 %
1521 % o map_id: ID of the map to look for in XML list.
1522 %
1523 % o exception: return any errors or warnings in this structure.
1524 %
1525 */
1526 static ThresholdMap *GetThresholdMapFile(const char *xml,const char *filename,
1527  const char *map_id,ExceptionInfo *exception)
1528 {
1529  char
1530  *p;
1531 
1532  const char
1533  *attribute,
1534  *content;
1535 
1536  double
1537  value;
1538 
1539  ssize_t
1540  i;
1541 
1542  ThresholdMap
1543  *map;
1544 
1545  XMLTreeInfo
1546  *description,
1547  *levels,
1548  *threshold,
1549  *thresholds;
1550 
1551  (void) LogMagickEvent(ConfigureEvent,GetMagickModule(),
1552  "Loading threshold map file \"%s\" ...",filename);
1553  map=(ThresholdMap *) NULL;
1554  thresholds=NewXMLTree(xml,exception);
1555  if (thresholds == (XMLTreeInfo *) NULL)
1556  return(map);
1557  for (threshold=GetXMLTreeChild(thresholds,"threshold");
1558  threshold != (XMLTreeInfo *) NULL;
1559  threshold=GetNextXMLTreeTag(threshold))
1560  {
1561  attribute=GetXMLTreeAttribute(threshold,"map");
1562  if ((attribute != (char *) NULL) && (LocaleCompare(map_id,attribute) == 0))
1563  break;
1564  attribute=GetXMLTreeAttribute(threshold,"alias");
1565  if ((attribute != (char *) NULL) && (LocaleCompare(map_id,attribute) == 0))
1566  break;
1567  }
1568  if (threshold == (XMLTreeInfo *) NULL)
1569  {
1570  thresholds=DestroyXMLTree(thresholds);
1571  return(map);
1572  }
1573  description=GetXMLTreeChild(threshold,"description");
1574  if (description == (XMLTreeInfo *) NULL)
1575  {
1576  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1577  "XmlMissingElement", "<description>, map \"%s\"",map_id);
1578  thresholds=DestroyXMLTree(thresholds);
1579  return(map);
1580  }
1581  levels=GetXMLTreeChild(threshold,"levels");
1582  if (levels == (XMLTreeInfo *) NULL)
1583  {
1584  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1585  "XmlMissingElement", "<levels>, map \"%s\"", map_id);
1586  thresholds=DestroyXMLTree(thresholds);
1587  return(map);
1588  }
1589  map=(ThresholdMap *) AcquireCriticalMemory(sizeof(*map));
1590  map->map_id=(char *) NULL;
1591  map->description=(char *) NULL;
1592  map->levels=(ssize_t *) NULL;
1593  attribute=GetXMLTreeAttribute(threshold,"map");
1594  if (attribute != (char *) NULL)
1595  map->map_id=ConstantString(attribute);
1596  content=GetXMLTreeContent(description);
1597  if (content != (char *) NULL)
1598  map->description=ConstantString(content);
1599  attribute=GetXMLTreeAttribute(levels,"width");
1600  if (attribute == (char *) NULL)
1601  {
1602  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1603  "XmlMissingAttribute", "<levels width>, map \"%s\"",map_id);
1604  thresholds=DestroyXMLTree(thresholds);
1605  map=DestroyThresholdMap(map);
1606  return(map);
1607  }
1608  map->width=StringToUnsignedLong(attribute);
1609  if (map->width == 0)
1610  {
1611  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1612  "XmlInvalidAttribute", "<levels width>, map \"%s\"",map_id);
1613  thresholds=DestroyXMLTree(thresholds);
1614  map=DestroyThresholdMap(map);
1615  return(map);
1616  }
1617  attribute=GetXMLTreeAttribute(levels,"height");
1618  if (attribute == (char *) NULL)
1619  {
1620  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1621  "XmlMissingAttribute", "<levels height>, map \"%s\"",map_id);
1622  thresholds=DestroyXMLTree(thresholds);
1623  map=DestroyThresholdMap(map);
1624  return(map);
1625  }
1626  map->height=StringToUnsignedLong(attribute);
1627  if (map->height == 0)
1628  {
1629  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1630  "XmlInvalidAttribute", "<levels height>, map \"%s\"",map_id);
1631  thresholds=DestroyXMLTree(thresholds);
1632  map=DestroyThresholdMap(map);
1633  return(map);
1634  }
1635  attribute=GetXMLTreeAttribute(levels,"divisor");
1636  if (attribute == (char *) NULL)
1637  {
1638  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1639  "XmlMissingAttribute", "<levels divisor>, map \"%s\"",map_id);
1640  thresholds=DestroyXMLTree(thresholds);
1641  map=DestroyThresholdMap(map);
1642  return(map);
1643  }
1644  map->divisor=(ssize_t) StringToLong(attribute);
1645  if (map->divisor < 2)
1646  {
1647  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1648  "XmlInvalidAttribute", "<levels divisor>, map \"%s\"",map_id);
1649  thresholds=DestroyXMLTree(thresholds);
1650  map=DestroyThresholdMap(map);
1651  return(map);
1652  }
1653  content=GetXMLTreeContent(levels);
1654  if (content == (char *) NULL)
1655  {
1656  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1657  "XmlMissingContent", "<levels>, map \"%s\"",map_id);
1658  thresholds=DestroyXMLTree(thresholds);
1659  map=DestroyThresholdMap(map);
1660  return(map);
1661  }
1662  map->levels=(ssize_t *) AcquireQuantumMemory((size_t) map->width,map->height*
1663  sizeof(*map->levels));
1664  if (map->levels == (ssize_t *) NULL)
1665  ThrowFatalException(ResourceLimitFatalError,"UnableToAcquireThresholdMap");
1666  for (i=0; i < (ssize_t) (map->width*map->height); i++)
1667  {
1668  map->levels[i]=(ssize_t) strtol(content,&p,10);
1669  if (p == content)
1670  {
1671  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1672  "XmlInvalidContent", "<level> too few values, map \"%s\"",map_id);
1673  thresholds=DestroyXMLTree(thresholds);
1674  map=DestroyThresholdMap(map);
1675  return(map);
1676  }
1677  if ((map->levels[i] < 0) || (map->levels[i] > map->divisor))
1678  {
1679  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1680  "XmlInvalidContent", "<level> %.20g out of range, map \"%s\"",
1681  (double) map->levels[i],map_id);
1682  thresholds=DestroyXMLTree(thresholds);
1683  map=DestroyThresholdMap(map);
1684  return(map);
1685  }
1686  content=p;
1687  }
1688  value=(double) strtol(content,&p,10);
1689  (void) value;
1690  if (p != content)
1691  {
1692  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1693  "XmlInvalidContent", "<level> too many values, map \"%s\"",map_id);
1694  thresholds=DestroyXMLTree(thresholds);
1695  map=DestroyThresholdMap(map);
1696  return(map);
1697  }
1698  thresholds=DestroyXMLTree(thresholds);
1699  return(map);
1700 }
1701 
1702 /*
1703 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1704 % %
1705 % %
1706 % %
1707 + L i s t T h r e s h o l d M a p F i l e %
1708 % %
1709 % %
1710 % %
1711 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1712 %
1713 % ListThresholdMapFile() lists the threshold maps and their descriptions
1714 % in the given XML file data.
1715 %
1716 % The format of the ListThresholdMaps method is:
1717 %
1718 % MagickBooleanType ListThresholdMaps(FILE *file,const char*xml,
1719 % const char *filename,ExceptionInfo *exception)
1720 %
1721 % A description of each parameter follows.
1722 %
1723 % o file: An pointer to the output FILE.
1724 %
1725 % o xml: The threshold map list in XML format.
1726 %
1727 % o filename: The threshold map XML filename.
1728 %
1729 % o exception: return any errors or warnings in this structure.
1730 %
1731 */
1732 MagickBooleanType ListThresholdMapFile(FILE *file,const char *xml,
1733  const char *filename,ExceptionInfo *exception)
1734 {
1735  const char
1736  *alias,
1737  *content,
1738  *map;
1739 
1740  XMLTreeInfo
1741  *description,
1742  *threshold,
1743  *thresholds;
1744 
1745  assert( xml != (char *) NULL );
1746  assert( file != (FILE *) NULL );
1747  (void) LogMagickEvent(ConfigureEvent,GetMagickModule(),
1748  "Loading threshold map file \"%s\" ...",filename);
1749  thresholds=NewXMLTree(xml,exception);
1750  if ( thresholds == (XMLTreeInfo *) NULL )
1751  return(MagickFalse);
1752  (void) FormatLocaleFile(file,"%-16s %-12s %s\n","Map","Alias","Description");
1753  (void) FormatLocaleFile(file,
1754  "----------------------------------------------------\n");
1755  threshold=GetXMLTreeChild(thresholds,"threshold");
1756  for ( ; threshold != (XMLTreeInfo *) NULL;
1757  threshold=GetNextXMLTreeTag(threshold))
1758  {
1759  map=GetXMLTreeAttribute(threshold,"map");
1760  if (map == (char *) NULL)
1761  {
1762  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1763  "XmlMissingAttribute", "<map>");
1764  thresholds=DestroyXMLTree(thresholds);
1765  return(MagickFalse);
1766  }
1767  alias=GetXMLTreeAttribute(threshold,"alias");
1768  description=GetXMLTreeChild(threshold,"description");
1769  if (description == (XMLTreeInfo *) NULL)
1770  {
1771  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1772  "XmlMissingElement", "<description>, map \"%s\"",map);
1773  thresholds=DestroyXMLTree(thresholds);
1774  return(MagickFalse);
1775  }
1776  content=GetXMLTreeContent(description);
1777  if (content == (char *) NULL)
1778  {
1779  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1780  "XmlMissingContent", "<description>, map \"%s\"", map);
1781  thresholds=DestroyXMLTree(thresholds);
1782  return(MagickFalse);
1783  }
1784  (void) FormatLocaleFile(file,"%-16s %-12s %s\n",map,alias ? alias : "",
1785  content);
1786  }
1787  thresholds=DestroyXMLTree(thresholds);
1788  return(MagickTrue);
1789 }
1790 
1791 /*
1792 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1793 % %
1794 % %
1795 % %
1796 % L i s t T h r e s h o l d M a p s %
1797 % %
1798 % %
1799 % %
1800 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1801 %
1802 % ListThresholdMaps() lists the threshold maps and their descriptions
1803 % as defined by "threshold.xml" to a file.
1804 %
1805 % The format of the ListThresholdMaps method is:
1806 %
1807 % MagickBooleanType ListThresholdMaps(FILE *file,ExceptionInfo *exception)
1808 %
1809 % A description of each parameter follows.
1810 %
1811 % o file: An pointer to the output FILE.
1812 %
1813 % o exception: return any errors or warnings in this structure.
1814 %
1815 */
1816 MagickExport MagickBooleanType ListThresholdMaps(FILE *file,
1817  ExceptionInfo *exception)
1818 {
1819  const StringInfo
1820  *option;
1821 
1823  *options;
1824 
1825  MagickStatusType
1826  status;
1827 
1828  status=MagickTrue;
1829  if (file == (FILE *) NULL)
1830  file=stdout;
1831  options=GetConfigureOptions(ThresholdsFilename,exception);
1832  (void) FormatLocaleFile(file,
1833  "\n Threshold Maps for Ordered Dither Operations\n");
1834  option=(const StringInfo *) GetNextValueInLinkedList(options);
1835  while (option != (const StringInfo *) NULL)
1836  {
1837  (void) FormatLocaleFile(file,"\nPath: %s\n\n",GetStringInfoPath(option));
1838  status&=(MagickStatusType) ListThresholdMapFile(file,(const char *)
1839  GetStringInfoDatum(option),GetStringInfoPath(option),exception);
1840  option=(const StringInfo *) GetNextValueInLinkedList(options);
1841  }
1842  options=DestroyConfigureOptions(options);
1843  return(status != 0 ? MagickTrue : MagickFalse);
1844 }
1845 
1846 /*
1847 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1848 % %
1849 % %
1850 % %
1851 % O r d e r e d D i t h e r I m a g e %
1852 % %
1853 % %
1854 % %
1855 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1856 %
1857 % OrderedDitherImage() will perform a ordered dither based on a number
1858 % of pre-defined dithering threshold maps, but over multiple intensity
1859 % levels, which can be different for different channels, according to the
1860 % input argument.
1861 %
1862 % The format of the OrderedDitherImage method is:
1863 %
1864 % MagickBooleanType OrderedDitherImage(Image *image,
1865 % const char *threshold_map,ExceptionInfo *exception)
1866 %
1867 % A description of each parameter follows:
1868 %
1869 % o image: the image.
1870 %
1871 % o threshold_map: A string containing the name of the threshold dither
1872 % map to use, followed by zero or more numbers representing the number
1873 % of color levels to dither between.
1874 %
1875 % Any level number less than 2 will be equivalent to 2, and means only
1876 % binary dithering will be applied to each color channel.
1877 %
1878 % No numbers also means a 2 level (bitmap) dither will be applied to all
1879 % channels, while a single number is the number of levels applied to each
1880 % channel in sequence. More numbers will be applied in turn to each of
1881 % the color channels.
1882 %
1883 % For example: "o3x3,6" will generate a 6 level posterization of the
1884 % image with an ordered 3x3 diffused pixel dither being applied between
1885 % each level. While checker,8,8,4 will produce a 332 colormaped image
1886 % with only a single checkerboard hash pattern (50% grey) between each
1887 % color level, to basically double the number of color levels with
1888 % a bare minimum of dithering.
1889 %
1890 % o exception: return any errors or warnings in this structure.
1891 %
1892 */
1893 MagickExport MagickBooleanType OrderedDitherImage(Image *image,
1894  const char *threshold_map,ExceptionInfo *exception)
1895 {
1896 #define DitherImageTag "Dither/Image"
1897 
1898  CacheView
1899  *image_view;
1900 
1901  char
1902  token[MagickPathExtent];
1903 
1904  const char
1905  *p;
1906 
1907  double
1908  levels[CompositePixelChannel];
1909 
1910  MagickBooleanType
1911  status;
1912 
1913  MagickOffsetType
1914  progress;
1915 
1916  ssize_t
1917  i,
1918  y;
1919 
1920  ThresholdMap
1921  *map;
1922 
1923  assert(image != (Image *) NULL);
1924  assert(image->signature == MagickCoreSignature);
1925  assert(exception != (ExceptionInfo *) NULL);
1926  assert(exception->signature == MagickCoreSignature);
1927  if (IsEventLogging() != MagickFalse)
1928  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1929  if (threshold_map == (const char *) NULL)
1930  return(MagickTrue);
1931  p=(char *) threshold_map;
1932  while (((isspace((int) ((unsigned char) *p)) != 0) || (*p == ',')) &&
1933  (*p != '\0'))
1934  p++;
1935  threshold_map=p;
1936  while (((isspace((int) ((unsigned char) *p)) == 0) && (*p != ',')) &&
1937  (*p != '\0'))
1938  {
1939  if ((p-threshold_map) >= (MagickPathExtent-1))
1940  break;
1941  token[p-threshold_map]=(*p);
1942  p++;
1943  }
1944  token[p-threshold_map]='\0';
1945  map=GetThresholdMap(token,exception);
1946  if (map == (ThresholdMap *) NULL)
1947  {
1948  (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1949  "InvalidArgument","%s : '%s'","ordered-dither",threshold_map);
1950  return(MagickFalse);
1951  }
1952  for (i=0; i < MaxPixelChannels; i++)
1953  levels[i]=2.0;
1954  p=strchr((char *) threshold_map,',');
1955  if ((p != (char *) NULL) && (isdigit((int) ((unsigned char) *(++p))) != 0))
1956  {
1957  (void) GetNextToken(p,&p,MagickPathExtent,token);
1958  for (i=0; (i < MaxPixelChannels); i++)
1959  levels[i]=StringToDouble(token,(char **) NULL);
1960  for (i=0; (*p != '\0') && (i < MaxPixelChannels); i++)
1961  {
1962  (void) GetNextToken(p,&p,MagickPathExtent,token);
1963  if (*token == ',')
1964  (void) GetNextToken(p,&p,MagickPathExtent,token);
1965  levels[i]=StringToDouble(token,(char **) NULL);
1966  }
1967  }
1968  for (i=0; i < MaxPixelChannels; i++)
1969  if (fabs(levels[i]) >= 1)
1970  levels[i]-=1.0;
1971  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1972  return(MagickFalse);
1973  status=MagickTrue;
1974  progress=0;
1975  image_view=AcquireAuthenticCacheView(image,exception);
1976 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1977  #pragma omp parallel for schedule(static) shared(progress,status) \
1978  magick_number_threads(image,image,image->rows,2)
1979 #endif
1980  for (y=0; y < (ssize_t) image->rows; y++)
1981  {
1982  ssize_t
1983  x;
1984 
1985  Quantum
1986  *magick_restrict q;
1987 
1988  if (status == MagickFalse)
1989  continue;
1990  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1991  if (q == (Quantum *) NULL)
1992  {
1993  status=MagickFalse;
1994  continue;
1995  }
1996  for (x=0; x < (ssize_t) image->columns; x++)
1997  {
1998  ssize_t
1999  j,
2000  n;
2001 
2002  n=0;
2003  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
2004  {
2005  ssize_t
2006  level,
2007  threshold;
2008 
2009  PixelChannel channel = GetPixelChannelChannel(image,j);
2010  PixelTrait traits = GetPixelChannelTraits(image,channel);
2011  if ((traits & UpdatePixelTrait) == 0)
2012  continue;
2013  if (fabs(levels[n]) < MagickEpsilon)
2014  {
2015  n++;
2016  continue;
2017  }
2018  threshold=(ssize_t) (QuantumScale*(double) q[j]*(levels[n]*
2019  (map->divisor-1)+1));
2020  level=threshold/(map->divisor-1);
2021  threshold-=level*(map->divisor-1);
2022  q[j]=ClampToQuantum((double) (level+(threshold >=
2023  map->levels[(x % (ssize_t) map->width)+(ssize_t) map->width*
2024  (y % (ssize_t) map->height)]))*(double) QuantumRange/levels[n]);
2025  n++;
2026  }
2027  q+=(ptrdiff_t) GetPixelChannels(image);
2028  }
2029  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2030  status=MagickFalse;
2031  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2032  {
2033  MagickBooleanType
2034  proceed;
2035 
2036 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2037  #pragma omp atomic
2038 #endif
2039  progress++;
2040  proceed=SetImageProgress(image,DitherImageTag,progress,image->rows);
2041  if (proceed == MagickFalse)
2042  status=MagickFalse;
2043  }
2044  }
2045  image_view=DestroyCacheView(image_view);
2046  map=DestroyThresholdMap(map);
2047  return(MagickTrue);
2048 }
2049 
2050 /*
2051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2052 % %
2053 % %
2054 % %
2055 % P e r c e p t i b l e I m a g e %
2056 % %
2057 % %
2058 % %
2059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2060 %
2061 % PerceptibleImage() set each pixel whose value is less than |epsilon| to
2062 % epsilon or -epsilon (whichever is closer) otherwise the pixel value remains
2063 % unchanged.
2064 %
2065 % The format of the PerceptibleImage method is:
2066 %
2067 % MagickBooleanType PerceptibleImage(Image *image,const double epsilon,
2068 % ExceptionInfo *exception)
2069 %
2070 % A description of each parameter follows:
2071 %
2072 % o image: the image.
2073 %
2074 % o epsilon: the epsilon threshold (e.g. 1.0e-9).
2075 %
2076 % o exception: return any errors or warnings in this structure.
2077 %
2078 */
2079 
2080 static inline Quantum PerceptibleThreshold(const Quantum quantum,
2081  const double epsilon)
2082 {
2083  double
2084  sign;
2085 
2086  sign=(double) quantum < 0.0 ? -1.0 : 1.0;
2087  if ((sign*(double) quantum) >= epsilon)
2088  return(quantum);
2089  return((Quantum) (sign*epsilon));
2090 }
2091 
2092 MagickExport MagickBooleanType PerceptibleImage(Image *image,
2093  const double epsilon,ExceptionInfo *exception)
2094 {
2095 #define PerceptibleImageTag "Perceptible/Image"
2096 
2097  CacheView
2098  *image_view;
2099 
2100  MagickBooleanType
2101  status;
2102 
2103  MagickOffsetType
2104  progress;
2105 
2106  ssize_t
2107  y;
2108 
2109  assert(image != (Image *) NULL);
2110  assert(image->signature == MagickCoreSignature);
2111  if (IsEventLogging() != MagickFalse)
2112  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2113  if (image->storage_class == PseudoClass)
2114  {
2115  ssize_t
2116  i;
2117 
2118  PixelInfo
2119  *magick_restrict q;
2120 
2121  q=image->colormap;
2122  for (i=0; i < (ssize_t) image->colors; i++)
2123  {
2124  if ((GetPixelChannelTraits(image,RedPixelChannel) & UpdatePixelTrait) != 0)
2125  q->red=(MagickRealType) PerceptibleThreshold(ClampToQuantum(q->red),
2126  epsilon);
2127  if ((GetPixelChannelTraits(image,GreenPixelChannel) & UpdatePixelTrait) != 0)
2128  q->green=(MagickRealType) PerceptibleThreshold(ClampToQuantum(q->green),
2129  epsilon);
2130  if ((GetPixelChannelTraits(image,BluePixelChannel) & UpdatePixelTrait) != 0)
2131  q->blue=(MagickRealType) PerceptibleThreshold(ClampToQuantum(q->blue),
2132  epsilon);
2133  if ((GetPixelChannelTraits(image,AlphaPixelChannel) & UpdatePixelTrait) != 0)
2134  q->alpha=(MagickRealType) PerceptibleThreshold(ClampToQuantum(q->alpha),
2135  epsilon);
2136  q++;
2137  }
2138  return(SyncImage(image,exception));
2139  }
2140  /*
2141  Perceptible image.
2142  */
2143  status=MagickTrue;
2144  progress=0;
2145  image_view=AcquireAuthenticCacheView(image,exception);
2146 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2147  #pragma omp parallel for schedule(static) shared(progress,status) \
2148  magick_number_threads(image,image,image->rows,2)
2149 #endif
2150  for (y=0; y < (ssize_t) image->rows; y++)
2151  {
2152  ssize_t
2153  x;
2154 
2155  Quantum
2156  *magick_restrict q;
2157 
2158  if (status == MagickFalse)
2159  continue;
2160  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2161  if (q == (Quantum *) NULL)
2162  {
2163  status=MagickFalse;
2164  continue;
2165  }
2166  for (x=0; x < (ssize_t) image->columns; x++)
2167  {
2168  ssize_t
2169  i;
2170 
2171  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2172  {
2173  PixelChannel channel = GetPixelChannelChannel(image,i);
2174  PixelTrait traits = GetPixelChannelTraits(image,channel);
2175  if ((traits & UpdatePixelTrait) != 0)
2176  q[i]=PerceptibleThreshold(q[i],epsilon);
2177  }
2178  q+=(ptrdiff_t) GetPixelChannels(image);
2179  }
2180  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2181  status=MagickFalse;
2182  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2183  {
2184  MagickBooleanType
2185  proceed;
2186 
2187 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2188  #pragma omp atomic
2189 #endif
2190  progress++;
2191  proceed=SetImageProgress(image,PerceptibleImageTag,progress,
2192  image->rows);
2193  if (proceed == MagickFalse)
2194  status=MagickFalse;
2195  }
2196  }
2197  image_view=DestroyCacheView(image_view);
2198  return(status);
2199 }
2200 
2201 /*
2202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2203 % %
2204 % %
2205 % %
2206 % R a n d o m T h r e s h o l d I m a g e %
2207 % %
2208 % %
2209 % %
2210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2211 %
2212 % RandomThresholdImage() changes the value of individual pixels based on the
2213 % intensity of each pixel compared to a random threshold. The result is a
2214 % low-contrast, two color image.
2215 %
2216 % The format of the RandomThresholdImage method is:
2217 %
2218 % MagickBooleanType RandomThresholdImage(Image *image,
2219 % const char *thresholds,ExceptionInfo *exception)
2220 %
2221 % A description of each parameter follows:
2222 %
2223 % o image: the image.
2224 %
2225 % o low,high: Specify the high and low thresholds. These values range from
2226 % 0 to QuantumRange.
2227 %
2228 % o exception: return any errors or warnings in this structure.
2229 %
2230 */
2231 MagickExport MagickBooleanType RandomThresholdImage(Image *image,
2232  const double min_threshold, const double max_threshold,ExceptionInfo *exception)
2233 {
2234 #define ThresholdImageTag "Threshold/Image"
2235 
2236  CacheView
2237  *image_view;
2238 
2239  MagickBooleanType
2240  status;
2241 
2242  MagickOffsetType
2243  progress;
2244 
2245  RandomInfo
2246  **magick_restrict random_info;
2247 
2248  ssize_t
2249  y;
2250 
2251 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2252  unsigned long
2253  key;
2254 #endif
2255 
2256  assert(image != (Image *) NULL);
2257  assert(image->signature == MagickCoreSignature);
2258  assert(exception != (ExceptionInfo *) NULL);
2259  assert(exception->signature == MagickCoreSignature);
2260  if (IsEventLogging() != MagickFalse)
2261  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2262  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2263  return(MagickFalse);
2264  /*
2265  Random threshold image.
2266  */
2267  status=MagickTrue;
2268  progress=0;
2269  random_info=AcquireRandomInfoTLS();
2270  image_view=AcquireAuthenticCacheView(image,exception);
2271 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2272  key=GetRandomSecretKey(random_info[0]);
2273  #pragma omp parallel for schedule(static) shared(progress,status) \
2274  magick_number_threads(image,image,image->rows,key == ~0UL ? 0 : 1)
2275 #endif
2276  for (y=0; y < (ssize_t) image->rows; y++)
2277  {
2278  const int
2279  id = GetOpenMPThreadId();
2280 
2281  Quantum
2282  *magick_restrict q;
2283 
2284  ssize_t
2285  x;
2286 
2287  if (status == MagickFalse)
2288  continue;
2289  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2290  if (q == (Quantum *) NULL)
2291  {
2292  status=MagickFalse;
2293  continue;
2294  }
2295  for (x=0; x < (ssize_t) image->columns; x++)
2296  {
2297  ssize_t
2298  i;
2299 
2300  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2301  {
2302  double
2303  threshold;
2304 
2305  PixelChannel channel = GetPixelChannelChannel(image,i);
2306  PixelTrait traits = GetPixelChannelTraits(image,channel);
2307  if ((traits & UpdatePixelTrait) == 0)
2308  continue;
2309  if ((double) q[i] < min_threshold)
2310  threshold=min_threshold;
2311  else
2312  if ((double) q[i] > max_threshold)
2313  threshold=max_threshold;
2314  else
2315  threshold=((double) QuantumRange*
2316  GetPseudoRandomValue(random_info[id]));
2317  q[i]=(double) q[i] <= threshold ? 0 : QuantumRange;
2318  }
2319  q+=(ptrdiff_t) GetPixelChannels(image);
2320  }
2321  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2322  status=MagickFalse;
2323  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2324  {
2325  MagickBooleanType
2326  proceed;
2327 
2328 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2329  #pragma omp atomic
2330 #endif
2331  progress++;
2332  proceed=SetImageProgress(image,ThresholdImageTag,progress,
2333  image->rows);
2334  if (proceed == MagickFalse)
2335  status=MagickFalse;
2336  }
2337  }
2338  image_view=DestroyCacheView(image_view);
2339  random_info=DestroyRandomInfoTLS(random_info);
2340  return(status);
2341 }
2342 
2343 /*
2344 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2345 % %
2346 % %
2347 % %
2348 % R a n g e T h r e s h o l d I m a g e %
2349 % %
2350 % %
2351 % %
2352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2353 %
2354 % RangeThresholdImage() applies soft and hard thresholding.
2355 %
2356 % The format of the RangeThresholdImage method is:
2357 %
2358 % MagickBooleanType RangeThresholdImage(Image *image,
2359 % const double low_black,const double low_white,const double high_white,
2360 % const double high_black,ExceptionInfo *exception)
2361 %
2362 % A description of each parameter follows:
2363 %
2364 % o image: the image.
2365 %
2366 % o low_black: Define the minimum black threshold value.
2367 %
2368 % o low_white: Define the minimum white threshold value.
2369 %
2370 % o high_white: Define the maximum white threshold value.
2371 %
2372 % o high_black: Define the maximum black threshold value.
2373 %
2374 % o exception: return any errors or warnings in this structure.
2375 %
2376 */
2377 MagickExport MagickBooleanType RangeThresholdImage(Image *image,
2378  const double low_black,const double low_white,const double high_white,
2379  const double high_black,ExceptionInfo *exception)
2380 {
2381 #define ThresholdImageTag "Threshold/Image"
2382 
2383  CacheView
2384  *image_view;
2385 
2386  MagickBooleanType
2387  status;
2388 
2389  MagickOffsetType
2390  progress;
2391 
2392  ssize_t
2393  y;
2394 
2395  assert(image != (Image *) NULL);
2396  assert(image->signature == MagickCoreSignature);
2397  if (IsEventLogging() != MagickFalse)
2398  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2399  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2400  return(MagickFalse);
2401  if (IsGrayColorspace(image->colorspace) != MagickFalse)
2402  (void) TransformImageColorspace(image,sRGBColorspace,exception);
2403  /*
2404  Range threshold image.
2405  */
2406  status=MagickTrue;
2407  progress=0;
2408  image_view=AcquireAuthenticCacheView(image,exception);
2409 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2410  #pragma omp parallel for schedule(static) shared(progress,status) \
2411  magick_number_threads(image,image,image->rows,2)
2412 #endif
2413  for (y=0; y < (ssize_t) image->rows; y++)
2414  {
2415  ssize_t
2416  x;
2417 
2418  Quantum
2419  *magick_restrict q;
2420 
2421  if (status == MagickFalse)
2422  continue;
2423  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2424  if (q == (Quantum *) NULL)
2425  {
2426  status=MagickFalse;
2427  continue;
2428  }
2429  for (x=0; x < (ssize_t) image->columns; x++)
2430  {
2431  double
2432  pixel;
2433 
2434  ssize_t
2435  i;
2436 
2437  pixel=GetPixelIntensity(image,q);
2438  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2439  {
2440  PixelChannel channel = GetPixelChannelChannel(image,i);
2441  PixelTrait traits = GetPixelChannelTraits(image,channel);
2442  if ((traits & UpdatePixelTrait) == 0)
2443  continue;
2444  if (image->channel_mask != AllChannels)
2445  pixel=(double) q[i];
2446  if (pixel < low_black)
2447  q[i]=(Quantum) 0;
2448  else
2449  if ((pixel >= low_black) && (pixel < low_white))
2450  q[i]=ClampToQuantum((double) QuantumRange*
2451  PerceptibleReciprocal(low_white-low_black)*(pixel-low_black));
2452  else
2453  if ((pixel >= low_white) && (pixel <= high_white))
2454  q[i]=QuantumRange;
2455  else
2456  if ((pixel > high_white) && (pixel <= high_black))
2457  q[i]=ClampToQuantum((double) QuantumRange*(double)
2458  PerceptibleReciprocal(high_black-high_white)*
2459  (high_black-pixel));
2460  else
2461  if (pixel > high_black)
2462  q[i]=(Quantum) 0;
2463  else
2464  q[i]=(Quantum) 0;
2465  }
2466  q+=(ptrdiff_t) GetPixelChannels(image);
2467  }
2468  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2469  status=MagickFalse;
2470  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2471  {
2472  MagickBooleanType
2473  proceed;
2474 
2475 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2476  #pragma omp atomic
2477 #endif
2478  progress++;
2479  proceed=SetImageProgress(image,ThresholdImageTag,progress,
2480  image->rows);
2481  if (proceed == MagickFalse)
2482  status=MagickFalse;
2483  }
2484  }
2485  image_view=DestroyCacheView(image_view);
2486  return(status);
2487 }
2488 
2489 /*
2490 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2491 % %
2492 % %
2493 % %
2494 % W h i t e T h r e s h o l d I m a g e %
2495 % %
2496 % %
2497 % %
2498 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2499 %
2500 % WhiteThresholdImage() is like ThresholdImage() but forces all pixels above
2501 % the threshold into white while leaving all pixels at or below the threshold
2502 % unchanged.
2503 %
2504 % The format of the WhiteThresholdImage method is:
2505 %
2506 % MagickBooleanType WhiteThresholdImage(Image *image,
2507 % const char *threshold,ExceptionInfo *exception)
2508 %
2509 % A description of each parameter follows:
2510 %
2511 % o image: the image.
2512 %
2513 % o threshold: Define the threshold value.
2514 %
2515 % o exception: return any errors or warnings in this structure.
2516 %
2517 */
2518 MagickExport MagickBooleanType WhiteThresholdImage(Image *image,
2519  const char *thresholds,ExceptionInfo *exception)
2520 {
2521 #define ThresholdImageTag "Threshold/Image"
2522 
2523  CacheView
2524  *image_view;
2525 
2526  GeometryInfo
2527  geometry_info;
2528 
2529  MagickBooleanType
2530  status;
2531 
2532  MagickOffsetType
2533  progress;
2534 
2535  PixelInfo
2536  threshold;
2537 
2538  MagickStatusType
2539  flags;
2540 
2541  ssize_t
2542  y;
2543 
2544  assert(image != (Image *) NULL);
2545  assert(image->signature == MagickCoreSignature);
2546  if (IsEventLogging() != MagickFalse)
2547  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2548  if (thresholds == (const char *) NULL)
2549  return(MagickTrue);
2550  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2551  return(MagickFalse);
2552  if (IsGrayColorspace(image->colorspace) != MagickFalse)
2553  (void) TransformImageColorspace(image,sRGBColorspace,exception);
2554  GetPixelInfo(image,&threshold);
2555  flags=ParseGeometry(thresholds,&geometry_info);
2556  threshold.red=geometry_info.rho;
2557  threshold.green=geometry_info.rho;
2558  threshold.blue=geometry_info.rho;
2559  threshold.black=geometry_info.rho;
2560  threshold.alpha=100.0;
2561  if ((flags & SigmaValue) != 0)
2562  threshold.green=geometry_info.sigma;
2563  if ((flags & XiValue) != 0)
2564  threshold.blue=geometry_info.xi;
2565  if ((flags & PsiValue) != 0)
2566  threshold.alpha=geometry_info.psi;
2567  if (threshold.colorspace == CMYKColorspace)
2568  {
2569  if ((flags & PsiValue) != 0)
2570  threshold.black=geometry_info.psi;
2571  if ((flags & ChiValue) != 0)
2572  threshold.alpha=geometry_info.chi;
2573  }
2574  if ((flags & PercentValue) != 0)
2575  {
2576  threshold.red*=((MagickRealType) QuantumRange/100.0);
2577  threshold.green*=((MagickRealType) QuantumRange/100.0);
2578  threshold.blue*=((MagickRealType) QuantumRange/100.0);
2579  threshold.black*=((MagickRealType) QuantumRange/100.0);
2580  threshold.alpha*=((MagickRealType) QuantumRange/100.0);
2581  }
2582  /*
2583  White threshold image.
2584  */
2585  status=MagickTrue;
2586  progress=0;
2587  image_view=AcquireAuthenticCacheView(image,exception);
2588 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2589  #pragma omp parallel for schedule(static) shared(progress,status) \
2590  magick_number_threads(image,image,image->rows,2)
2591 #endif
2592  for (y=0; y < (ssize_t) image->rows; y++)
2593  {
2594  ssize_t
2595  x;
2596 
2597  Quantum
2598  *magick_restrict q;
2599 
2600  if (status == MagickFalse)
2601  continue;
2602  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2603  if (q == (Quantum *) NULL)
2604  {
2605  status=MagickFalse;
2606  continue;
2607  }
2608  for (x=0; x < (ssize_t) image->columns; x++)
2609  {
2610  double
2611  pixel;
2612 
2613  ssize_t
2614  i;
2615 
2616  pixel=GetPixelIntensity(image,q);
2617  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2618  {
2619  PixelChannel channel = GetPixelChannelChannel(image,i);
2620  PixelTrait traits = GetPixelChannelTraits(image,channel);
2621  if ((traits & UpdatePixelTrait) == 0)
2622  continue;
2623  if (image->channel_mask != AllChannels)
2624  pixel=(double) q[i];
2625  if (pixel > GetPixelInfoChannel(&threshold,channel))
2626  q[i]=QuantumRange;
2627  }
2628  q+=(ptrdiff_t) GetPixelChannels(image);
2629  }
2630  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2631  status=MagickFalse;
2632  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2633  {
2634  MagickBooleanType
2635  proceed;
2636 
2637 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2638  #pragma omp atomic
2639 #endif
2640  progress++;
2641  proceed=SetImageProgress(image,ThresholdImageTag,progress,image->rows);
2642  if (proceed == MagickFalse)
2643  status=MagickFalse;
2644  }
2645  }
2646  image_view=DestroyCacheView(image_view);
2647  return(status);
2648 }
_GeometryInfo
Definition: geometry.h:105
_CacheView
Definition: cache-view.c:65
_XMLTreeInfo
Definition: xml-tree.c:77
_Image
Definition: image.h:131
_PixelInfo
Definition: pixel.h:181
_LinkedListInfo
Definition: linked-list.c:60
_ExceptionInfo
Definition: exception.h:101
_ThresholdMap
Definition: threshold.c:101
_StringInfo
Definition: string_.h:27
_RandomInfo
Definition: random.c:83