← Back to branch summary

~ubuntu-branches/ubuntu/quantal/imagemagick/quantal

« back to all changes in this revision

Viewing changes to magick/statistic.c

  • Committer: Bazaar Package Importer
  • Author(s): Muharem Hrnjadovic
  • Date: 2009-06-04 13:01:13 UTC
  • mfrom: (1.1.5 upstream) (6.1.1 sid)
  • Revision ID: james.westby@ubuntu.com-20090604130113-my9114jxmafpwew3
Tags: 7:6.5.1.0-1.1ubuntu1
* Merge from debian unstable, remaining changes:
  - (Build-)depend on libltdl7-dev instead of libltdl3-dev (the armel buildds
    currently have both available).
  - Don't build-dep on librsvg, it brings in excessive dependencies

Show diffs side-by-side

added added

removed removed

Lines of Context:
magick/statistic.c
17
17
%                                 July 1992                                   %
18
18
%                                                                             %
19
19
%                                                                             %
20
 
%  Copyright 1999-2008 ImageMagick Studio LLC, a non-profit organization      %
 
20
%  Copyright 1999-2009 ImageMagick Studio LLC, a non-profit organization      %
21
21
%  dedicated to making software imaging solutions freely available.           %
22
22
%                                                                             %
23
23
%  You may not use this file except in compliance with the License.  You may  %
112
112
%
113
113
%    o image: the image.
114
114
%
115
 
%    o exception: Return any errors or warnings in this structure.
 
115
%    o exception: return any errors or warnings in this structure.
116
116
%
117
117
*/
118
118
MagickExport RectangleInfo GetImageBoundingBox(const Image *image,
135
135
    *p;
136
136
 
137
137
  ViewInfo
138
 
    **image_view;
 
138
    *image_view;
139
139
 
140
140
  assert(image != (Image *) NULL);
141
141
  assert(image->signature == MagickSignature);
146
146
  bounds.x=(long) image->columns;
147
147
  bounds.y=(long) image->rows;
148
148
  GetMagickPixelPacket(image,&target[0]);
149
 
  image_view=AcquireCacheViewThreadSet(image);
150
 
  p=AcquireCacheViewPixels(image_view[0],0,0,1,1,exception);
 
149
  image_view=AcquireCacheView(image);
 
150
  p=GetCacheViewVirtualPixels(image_view,0,0,1,1,exception);
151
151
  if (p == (const PixelPacket *) NULL)
152
152
    {
153
 
      image_view=DestroyCacheViewThreadSet(image_view);
 
153
      image_view=DestroyCacheView(image_view);
154
154
      return(bounds);
155
155
    }
156
 
  SetMagickPixelPacket(image,p,GetCacheViewIndexes(image_view[0]),&target[0]);
 
156
  SetMagickPixelPacket(image,p,GetCacheViewAuthenticIndexQueue(image_view),
 
157
    &target[0]);
157
158
  GetMagickPixelPacket(image,&target[1]);
158
 
  p=AcquireCacheViewPixels(image_view[0],(long) image->columns-1,0,1,1,
 
159
  p=GetCacheViewVirtualPixels(image_view,(long) image->columns-1,0,1,1,
159
160
    exception);
160
 
  SetMagickPixelPacket(image,p,GetCacheViewIndexes(image_view[0]),&target[1]);
 
161
  SetMagickPixelPacket(image,p,GetCacheViewAuthenticIndexQueue(image_view),
 
162
    &target[1]);
161
163
  GetMagickPixelPacket(image,&target[2]);
162
 
  p=AcquireCacheViewPixels(image_view[0],0,(long) image->rows-1,1,1,exception);
163
 
  SetMagickPixelPacket(image,p,GetCacheViewIndexes(image_view[0]),&target[2]);
 
164
  p=GetCacheViewVirtualPixels(image_view,0,(long) image->rows-1,1,1,exception);
 
165
  SetMagickPixelPacket(image,p,GetCacheViewAuthenticIndexQueue(image_view),
 
166
    &target[2]);
164
167
  status=MagickTrue;
165
168
  GetMagickPixelPacket(image,&zero);
166
169
#if defined(MAGICKCORE_OPENMP_SUPPORT)
167
 
  #pragma omp parallel for schedule(static,64) shared(status)
 
170
  #pragma omp parallel for schedule(dynamic,4) shared(status)
168
171
#endif
169
172
  for (y=0; y < (long) image->rows; y++)
170
173
  {
171
174
    MagickPixelPacket
172
175
      pixel;
173
176
 
 
177
    RectangleInfo
 
178
      bounding_box;
 
179
 
 
180
    register const IndexPacket
 
181
      *indexes;
 
182
 
174
183
    register const PixelPacket
175
184
      *p;
176
185
 
177
 
    register IndexPacket
178
 
      *indexes;
179
 
 
180
186
    register long
181
 
      id,
182
187
      x;
183
188
 
184
189
    if (status == MagickFalse)
185
190
      continue;
186
 
    id=GetCacheViewThreadId();
187
 
    p=AcquireCacheViewPixels(image_view[id],0,y,image->columns,1,exception);
 
191
#if defined(HAVE_OPENMP)
 
192
#  pragma omp critical (MagickCore_GetImageBoundingBox)
 
193
#endif
 
194
    bounding_box=bounds;
 
195
    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
188
196
    if (p == (const PixelPacket *) NULL)
189
197
      {
190
198
        status=MagickFalse;
191
199
        continue;
192
200
      }
193
 
    indexes=GetCacheViewIndexes(image_view[id]);
 
201
    indexes=GetCacheViewVirtualIndexQueue(image_view);
194
202
    pixel=zero;
195
203
    for (x=0; x < (long) image->columns; x++)
196
204
    {
197
205
      SetMagickPixelPacket(image,p,indexes+x,&pixel);
198
 
      if ((x < bounds.x) &&
 
206
      if ((x < bounding_box.x) &&
199
207
          (IsMagickColorSimilar(&pixel,&target[0]) == MagickFalse))
200
 
        bounds.x=x;
201
 
      if ((x > (long) bounds.width) &&
 
208
        bounding_box.x=x;
 
209
      if ((x > (long) bounding_box.width) &&
202
210
          (IsMagickColorSimilar(&pixel,&target[1]) == MagickFalse))
203
 
        bounds.width=(unsigned long) x;
204
 
      if ((y < bounds.y) &&
 
211
        bounding_box.width=(unsigned long) x;
 
212
      if ((y < bounding_box.y) &&
205
213
          (IsMagickColorSimilar(&pixel,&target[0]) == MagickFalse))
206
 
        bounds.y=y;
207
 
      if ((y > (long) bounds.height) &&
 
214
        bounding_box.y=y;
 
215
      if ((y > (long) bounding_box.height) &&
208
216
          (IsMagickColorSimilar(&pixel,&target[2]) == MagickFalse))
209
 
        bounds.height=(unsigned long) y;
 
217
        bounding_box.height=(unsigned long) y;
210
218
      p++;
211
219
    }
 
220
#if defined(HAVE_OPENMP)
 
221
#  pragma omp critical (MagickCore_GetImageBoundingBox)
 
222
#endif
 
223
    {
 
224
      if (bounding_box.x < bounds.x)
 
225
        bounds.x=bounding_box.x;
 
226
      if (bounding_box.y < bounds.y)
 
227
        bounds.y=bounding_box.y;
 
228
      if (bounding_box.width > bounds.width)
 
229
        bounds.width=bounding_box.width;
 
230
      if (bounding_box.height > bounds.height)
 
231
        bounds.height=bounding_box.height;
 
232
    }
212
233
  }
213
 
  image_view=DestroyCacheViewThreadSet(image_view);
 
234
  image_view=DestroyCacheView(image_view);
214
235
  if ((bounds.width == 0) || (bounds.height == 0))
215
236
    (void) ThrowMagickException(exception,GetMagickModule(),OptionWarning,
216
237
      "GeometryDoesNotContainImage","`%s'",image->filename);
247
268
%
248
269
%    o channel: the channel.
249
270
%
250
 
%    o exception: Return any errors or warnings in this structure.
 
271
%    o exception: return any errors or warnings in this structure.
251
272
%
252
273
*/
253
274
 
271
292
 
272
293
  unsigned long
273
294
    *current_depth,
274
 
    depth;
275
 
    
 
295
    depth,
 
296
    number_threads;
 
297
 
276
298
  ViewInfo
277
 
    **image_view;
278
 
  
 
299
    *image_view;
 
300
 
279
301
  /*
280
302
    Compute image depth.
281
303
  */
283
305
  assert(image->signature == MagickSignature);
284
306
  if (image->debug != MagickFalse)
285
307
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
286
 
  current_depth=(unsigned long *) AcquireQuantumMemory(
287
 
    GetCacheViewMaximumThreads(),sizeof(*current_depth));
 
308
  number_threads=GetPixelCacheMaximumThreads();
 
309
  current_depth=(unsigned long *) AcquireQuantumMemory(number_threads,
 
310
    sizeof(*current_depth));
288
311
  if (current_depth == (unsigned long *) NULL)
289
312
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
290
313
  status=MagickTrue;
291
 
  for (id=0; id < (long) GetCacheViewMaximumThreads(); id++)
 
314
  for (id=0; id < (long) number_threads; id++)
292
315
    current_depth[id]=1;
293
316
  if ((image->storage_class == PseudoClass) && (image->matte == MagickFalse))
294
317
    {
300
323
 
301
324
      p=image->colormap;
302
325
#if defined(MAGICKCORE_OPENMP_SUPPORT)
303
 
      #pragma omp parallel for schedule(static,64) shared(status)
 
326
  #pragma omp parallel for schedule(dynamic,4) shared(status)
304
327
#endif
305
328
      for (i=0; i < (long) image->colors; i++)
306
 
      {  
 
329
      {
307
330
        if (status == MagickFalse)
308
331
          continue;
309
 
        id=GetCacheViewThreadId();
 
332
        id=GetPixelCacheThreadId();
310
333
        while (current_depth[id] < MAGICKCORE_QUANTUM_DEPTH)
311
334
        {
312
335
          MagickStatusType
313
336
            status;
314
337
 
315
338
          QuantumAny
316
 
            scale;
 
339
            range;
317
340
 
318
341
          status=0;
319
 
          scale=GetQuantumScale(current_depth[id]);
 
342
          range=GetQuantumRange(current_depth[id]);
320
343
          if ((channel & RedChannel) != 0)
321
344
            status|=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,
322
 
              current_depth[id],scale),current_depth[id],scale);
 
345
              range),range);
323
346
          if ((channel & GreenChannel) != 0)
324
347
            status|=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
325
 
              current_depth[id],scale),current_depth[id],scale);
 
348
              range),range);
326
349
          if ((channel & BlueChannel) != 0)
327
350
            status|=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
328
 
              current_depth[id],scale),current_depth[id],scale);
 
351
              range),range);
329
352
          if (status == 0)
330
353
            break;
331
354
          current_depth[id]++;
333
356
        p++;
334
357
      }
335
358
      depth=current_depth[0];
336
 
      for (id=1; id < (long) GetCacheViewMaximumThreads(); id++)
 
359
      for (id=1; id < (long) number_threads; id++)
337
360
        if (depth < current_depth[id])
338
361
          depth=current_depth[id];
339
362
      current_depth=(unsigned long *) RelinquishMagickMemory(current_depth);
340
363
      return(depth);
341
364
    }
342
 
  image_view=AcquireCacheViewThreadSet(image);
 
365
  image_view=AcquireCacheView(image);
343
366
#if defined(MAGICKCORE_OPENMP_SUPPORT)
344
 
  #pragma omp parallel for schedule(static,64) shared(status)
 
367
  #pragma omp parallel for schedule(dynamic,4) shared(status)
345
368
#endif
346
369
  for (y=0; y < (long) image->rows; y++)
347
370
  {
357
380
 
358
381
    if (status == MagickFalse)
359
382
      continue;
360
 
    id=GetCacheViewThreadId();
361
 
    p=AcquireCacheViewPixels(image_view[id],0,y,image->columns,1,exception);
 
383
    id=GetPixelCacheThreadId();
 
384
    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
362
385
    if (p == (const PixelPacket *) NULL)
363
386
      continue;
364
 
    indexes=AcquireCacheViewIndexes(image_view[id]);
 
387
    indexes=GetCacheViewVirtualIndexQueue(image_view);
365
388
    for (x=0; x < (long) image->columns; x++)
366
 
    {  
 
389
    {
367
390
      while (current_depth[id] < MAGICKCORE_QUANTUM_DEPTH)
368
391
      {
369
392
        MagickStatusType
370
393
          status;
371
394
 
372
395
        QuantumAny
373
 
          scale;
 
396
          range;
374
397
 
375
398
        status=0;
376
 
        scale=GetQuantumScale(current_depth[id]);
 
399
        range=GetQuantumRange(current_depth[id]);
377
400
        if ((channel & RedChannel) != 0)
378
 
          status|=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,
379
 
            current_depth[id],scale),current_depth[id],scale);
 
401
          status|=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),
 
402
            range);
380
403
        if ((channel & GreenChannel) != 0)
381
404
          status|=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
382
 
            current_depth[id],scale),current_depth[id],scale);
 
405
            range),range);
383
406
        if ((channel & BlueChannel) != 0)
384
 
          status|=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
385
 
            current_depth[id],scale),current_depth[id],scale);
 
407
          status|=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,range),
 
408
            range);
386
409
        if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
387
410
          status|=p->opacity != ScaleAnyToQuantum(ScaleQuantumToAny(p->opacity,
388
 
            current_depth[id],scale),current_depth[id],scale);
 
411
            range),range);
389
412
        if (((channel & IndexChannel) != 0) &&
390
413
            (image->colorspace == CMYKColorspace))
391
414
          status|=indexes[x] != ScaleAnyToQuantum(ScaleQuantumToAny(indexes[x],
392
 
            current_depth[id],scale),current_depth[id],scale);
 
415
            range),range);
393
416
        if (status == 0)
394
417
          break;
395
418
        current_depth[id]++;
399
422
    if (current_depth[id] == MAGICKCORE_QUANTUM_DEPTH)
400
423
      status=MagickFalse;
401
424
  }
402
 
  image_view=DestroyCacheViewThreadSet(image_view);
 
425
  image_view=DestroyCacheView(image_view);
403
426
  depth=current_depth[0];
404
 
  for (id=1; id < (long) GetCacheViewMaximumThreads(); id++)
 
427
  for (id=1; id < (long) number_threads; id++)
405
428
    if (depth < current_depth[id])
406
429
      depth=current_depth[id];
407
430
  current_depth=(unsigned long *) RelinquishMagickMemory(current_depth);
437
460
%
438
461
%    o maxima: the maximum value in the channel.
439
462
%
440
 
%    o exception: Return any errors or warnings in this structure.
 
463
%    o exception: return any errors or warnings in this structure.
441
464
%
442
465
*/
443
466
 
498
521
%
499
522
%    o standard_deviation: the standard deviation of the channel.
500
523
%
501
 
%    o exception: Return any errors or warnings in this structure.
 
524
%    o exception: return any errors or warnings in this structure.
502
525
%
503
526
*/
504
527
 
517
540
  const ChannelType channel,double *mean,double *standard_deviation,
518
541
  ExceptionInfo *exception)
519
542
{
520
 
#define PixelSquared(x)  ((x)*(x))
521
 
 
522
543
  double
523
544
    area;
524
545
 
543
564
  area=0.0;
544
565
  for (y=0; y < (long) image->rows; y++)
545
566
  {
546
 
    p=AcquireImagePixels(image,0,y,image->columns,1,exception);
 
567
    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
547
568
    if (p == (const PixelPacket *) NULL)
548
569
      break;
549
 
    indexes=AcquireIndexes(image);
 
570
    indexes=GetVirtualIndexQueue(image);
550
571
    for (x=0; x < (long) image->columns; x++)
551
572
    {
552
573
      if ((channel & RedChannel) != 0)
599
620
%                                                                             %
600
621
%                                                                             %
601
622
%                                                                             %
 
623
%   G e t I m a g e C h a n n e l K u r t o s i s                             %
 
624
%                                                                             %
 
625
%                                                                             %
 
626
%                                                                             %
 
627
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
628
%
 
629
%  GetImageChannelKurtosis() returns the kurtosis and skewness of one or more
 
630
%  image channels.
 
631
%
 
632
%  The format of the GetImageChannelKurtosis method is:
 
633
%
 
634
%      MagickBooleanType GetImageChannelKurtosis(const Image *image,
 
635
%        const ChannelType channel,double *kurtosis,double *skewness,
 
636
%        ExceptionInfo *exception)
 
637
%
 
638
%  A description of each parameter follows:
 
639
%
 
640
%    o image: the image.
 
641
%
 
642
%    o channel: the channel.
 
643
%
 
644
%    o kurtosis: the kurtosis of the channel.
 
645
%
 
646
%    o skewness: the skewness of the channel.
 
647
%
 
648
%    o exception: return any errors or warnings in this structure.
 
649
%
 
650
*/
 
651
 
 
652
MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
 
653
  double *kurtosis,double *skewness,ExceptionInfo *exception)
 
654
{
 
655
  MagickBooleanType
 
656
    status;
 
657
 
 
658
  status=GetImageChannelKurtosis(image,AllChannels,kurtosis,skewness,
 
659
    exception);
 
660
  return(status);
 
661
}
 
662
 
 
663
MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
 
664
  const ChannelType channel,double *kurtosis,double *skewness,
 
665
  ExceptionInfo *exception)
 
666
{
 
667
  double
 
668
    area,
 
669
    mean,
 
670
    standard_deviation,
 
671
    sum_squares,
 
672
    sum_cubes,
 
673
    sum_fourth_power;
 
674
 
 
675
  long
 
676
    y;
 
677
 
 
678
  register const IndexPacket
 
679
    *indexes;
 
680
 
 
681
  register const PixelPacket
 
682
    *p;
 
683
 
 
684
  register long
 
685
    x;
 
686
 
 
687
  assert(image != (Image *) NULL);
 
688
  assert(image->signature == MagickSignature);
 
689
  if (image->debug != MagickFalse)
 
690
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 
691
  *kurtosis=0.0;
 
692
  *skewness=0.0;
 
693
  area=0.0;
 
694
  mean=0.0;
 
695
  standard_deviation=0.0;
 
696
  sum_squares=0.0;
 
697
  sum_cubes=0.0;
 
698
  sum_fourth_power=0.0;
 
699
  for (y=0; y < (long) image->rows; y++)
 
700
  {
 
701
    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 
702
    if (p == (const PixelPacket *) NULL)
 
703
      break;
 
704
    indexes=GetVirtualIndexQueue(image);
 
705
    for (x=0; x < (long) image->columns; x++)
 
706
    {
 
707
      if ((channel & RedChannel) != 0)
 
708
        {
 
709
          mean+=p->red;
 
710
          sum_squares+=(double) p->red*p->red;
 
711
          sum_cubes+=(double) p->red*p->red*p->red;
 
712
          sum_fourth_power+=(double) p->red*p->red*p->red*p->red;
 
713
          area++;
 
714
        }
 
715
      if ((channel & GreenChannel) != 0)
 
716
        {
 
717
          mean+=p->green;
 
718
          sum_squares+=(double) p->green*p->green;
 
719
          sum_cubes+=(double) p->green*p->green*p->green;
 
720
          sum_fourth_power+=(double) p->green*p->green*p->green*p->green;
 
721
          area++;
 
722
        }
 
723
      if ((channel & BlueChannel) != 0)
 
724
        {
 
725
          mean+=p->blue;
 
726
          sum_squares+=(double) p->blue*p->blue;
 
727
          sum_cubes+=(double) p->blue*p->blue*p->blue;
 
728
          sum_fourth_power+=(double) p->blue*p->blue*p->blue*p->blue;
 
729
          area++;
 
730
        }
 
731
      if ((channel & OpacityChannel) != 0)
 
732
        {
 
733
          mean+=p->opacity;
 
734
          sum_squares+=(double) p->opacity*p->opacity;
 
735
          sum_cubes+=(double) p->opacity*p->opacity*p->opacity;
 
736
          sum_fourth_power+=(double) p->opacity*p->opacity*p->opacity*
 
737
            p->opacity;
 
738
          area++;
 
739
        }
 
740
      if (((channel & IndexChannel) != 0) &&
 
741
          (image->colorspace == CMYKColorspace))
 
742
        {
 
743
          mean+=indexes[x];
 
744
          sum_squares+=(double) indexes[x]*indexes[x];
 
745
          sum_cubes+=(double) indexes[x]*indexes[x]*indexes[x];
 
746
          sum_fourth_power+=(double) indexes[x]*indexes[x]*indexes[x]*
 
747
            indexes[x];
 
748
          area++;
 
749
        }
 
750
      p++;
 
751
    }
 
752
  }
 
753
  if (y < (long) image->rows)
 
754
    return(MagickFalse);
 
755
  if (area != 0.0)
 
756
    {
 
757
      mean/=area;
 
758
      sum_squares/=area;
 
759
      sum_cubes/=area;
 
760
      sum_fourth_power/=area;
 
761
    }
 
762
  standard_deviation=sqrt(sum_squares-(mean*mean));
 
763
  if (standard_deviation != 0.0)
 
764
    {
 
765
      *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
 
766
        3.0*mean*mean*mean*mean;
 
767
      *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
 
768
        standard_deviation;
 
769
      *kurtosis-=3.0;
 
770
      *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
 
771
      *skewness/=standard_deviation*standard_deviation*standard_deviation;
 
772
    }
 
773
  return(y == (long) image->rows ? MagickTrue : MagickFalse);
 
774
}
 
775
 
 
776
/*
 
777
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
778
%                                                                             %
 
779
%                                                                             %
 
780
%                                                                             %
602
781
%   G e t I m a g e C h a n n e l R a n g e                                   %
603
782
%                                                                             %
604
783
%                                                                             %
623
802
%
624
803
%    o maxima: the maximum value in the channel.
625
804
%
626
 
%    o exception: Return any errors or warnings in this structure.
 
805
%    o exception: return any errors or warnings in this structure.
627
806
%
628
807
*/
629
808
 
661
840
  GetMagickPixelPacket(image,&pixel);
662
841
  for (y=0; y < (long) image->rows; y++)
663
842
  {
664
 
    p=AcquireImagePixels(image,0,y,image->columns,1,exception);
 
843
    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
665
844
    if (p == (const PixelPacket *) NULL)
666
845
      break;
667
 
    indexes=AcquireIndexes(image);
 
846
    indexes=GetVirtualIndexQueue(image);
668
847
    for (x=0; x < (long) image->columns; x++)
669
848
    {
670
849
      SetMagickPixelPacket(image,p,indexes+x,&pixel);
722
901
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
723
902
%
724
903
%  GetImageChannelStatistics() returns statistics for each channel in the
725
 
%  image.  The statistics incude the channel depth, its minima, maxima, mean,
726
 
%  and the standard deviation.  You can access the red channel mean, for
727
 
%  example, like this:
 
904
%  image.  The statistics include the channel depth, its minima, maxima, mean,
 
905
%  standard deviation, kurtosis and skewness.  You can access the red channel
 
906
%  mean, for example, like this:
728
907
%
729
908
%      channel_statistics=GetImageChannelStatistics(image,excepton);
730
909
%      red_mean=channel_statistics[RedChannel].mean;
740
919
%
741
920
%    o image: the image.
742
921
%
743
 
%    o exception: Return any errors or warnings in this structure.
 
922
%    o exception: return any errors or warnings in this structure.
744
923
%
745
924
*/
746
925
 
765
944
    *channel_statistics;
766
945
 
767
946
  double
768
 
    area;
 
947
    area,
 
948
    sum_squares,
 
949
    sum_cubes;
769
950
 
770
951
  long
771
952
    y;
774
955
    status;
775
956
 
776
957
  QuantumAny
777
 
    scale;
 
958
    range;
778
959
 
779
960
  register const IndexPacket
780
961
    *indexes;
811
992
    channel_statistics[i].minima=1.0E+37;
812
993
    channel_statistics[i].mean=0.0;
813
994
    channel_statistics[i].standard_deviation=0.0;
 
995
    channel_statistics[i].kurtosis=0.0;
 
996
    channel_statistics[i].skewness=0.0;
814
997
  }
815
998
  y=(long) image->rows;
816
999
  for (y=0; y < (long) image->rows; y++)
817
1000
  {
818
 
    p=AcquireImagePixels(image,0,y,image->columns,1,exception);
 
1001
    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
819
1002
    if (p == (const PixelPacket *) NULL)
820
1003
      break;
821
 
    indexes=AcquireIndexes(image);
 
1004
    indexes=GetVirtualIndexQueue(image);
822
1005
    for (x=0; x < (long) image->columns; )
823
1006
    {
824
1007
      if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
825
1008
        {
826
1009
          depth=channel_statistics[RedChannel].depth;
827
 
          scale=GetQuantumScale(depth);
828
 
          status=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,depth,
829
 
            scale),depth,scale) ? MagickTrue : MagickFalse;
 
1010
          range=GetQuantumRange(depth);
 
1011
          status=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),
 
1012
            range) ? MagickTrue : MagickFalse;
830
1013
          if (status != MagickFalse)
831
1014
            {
832
1015
              channel_statistics[RedChannel].depth++;
836
1019
      if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
837
1020
        {
838
1021
          depth=channel_statistics[GreenChannel].depth;
839
 
          scale=GetQuantumScale(depth);
 
1022
          range=GetQuantumRange(depth);
840
1023
          status=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
841
 
            depth,scale),depth,scale) ? MagickTrue : MagickFalse;
 
1024
            range),range) ? MagickTrue : MagickFalse;
842
1025
          if (status != MagickFalse)
843
1026
            {
844
1027
              channel_statistics[GreenChannel].depth++;
848
1031
      if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
849
1032
        {
850
1033
          depth=channel_statistics[BlueChannel].depth;
851
 
          scale=GetQuantumScale(depth);
 
1034
          range=GetQuantumRange(depth);
852
1035
          status=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
853
 
            depth,scale),depth,scale) ? MagickTrue : MagickFalse;
 
1036
            range),range) ? MagickTrue : MagickFalse;
854
1037
          if (status != MagickFalse)
855
1038
            {
856
1039
              channel_statistics[BlueChannel].depth++;
860
1043
      if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
861
1044
        {
862
1045
          depth=channel_statistics[OpacityChannel].depth;
863
 
          scale=GetQuantumScale(depth);
 
1046
          range=GetQuantumRange(depth);
864
1047
          status=p->opacity != ScaleAnyToQuantum(ScaleQuantumToAny(p->opacity,
865
 
            depth,scale),depth,scale) ? MagickTrue : MagickFalse;
 
1048
            range),range) ? MagickTrue : MagickFalse;
866
1049
          if (status != MagickFalse)
867
1050
            {
868
1051
              channel_statistics[OpacityChannel].depth++;
874
1057
          if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
875
1058
            {
876
1059
              depth=channel_statistics[BlackChannel].depth;
877
 
              scale=GetQuantumScale(depth);
 
1060
              range=GetQuantumRange(depth);
878
1061
              status=indexes[x] != ScaleAnyToQuantum(ScaleQuantumToAny(
879
 
                indexes[x],depth,scale),depth,scale) ? MagickTrue : MagickFalse;
 
1062
                indexes[x],range),range) ? MagickTrue : MagickFalse;
880
1063
              if (status != MagickFalse)
881
1064
                {
882
1065
                  channel_statistics[BlackChannel].depth++;
889
1072
      if ((double) p->red > channel_statistics[RedChannel].maxima)
890
1073
        channel_statistics[RedChannel].maxima=(double) p->red;
891
1074
      channel_statistics[RedChannel].mean+=p->red;
892
 
      channel_statistics[RedChannel].standard_deviation+=(double)
 
1075
      channel_statistics[RedChannel].standard_deviation+=(double) p->red*p->red;
 
1076
      channel_statistics[RedChannel].kurtosis+=(double) p->red*p->red*
893
1077
        p->red*p->red;
 
1078
      channel_statistics[RedChannel].skewness+=(double) p->red*p->red*p->red;
894
1079
      if ((double) p->green < channel_statistics[GreenChannel].minima)
895
1080
        channel_statistics[GreenChannel].minima=(double) p->green;
896
1081
      if ((double) p->green > channel_statistics[GreenChannel].maxima)
897
1082
        channel_statistics[GreenChannel].maxima=(double) p->green;
898
1083
      channel_statistics[GreenChannel].mean+=p->green;
899
 
      channel_statistics[GreenChannel].standard_deviation+=(double)
 
1084
      channel_statistics[GreenChannel].standard_deviation+=(double) p->green*
 
1085
        p->green;
 
1086
      channel_statistics[GreenChannel].kurtosis+=(double) p->green*p->green*
900
1087
        p->green*p->green;
 
1088
      channel_statistics[GreenChannel].skewness+=(double) p->green*p->green*
 
1089
        p->green;
901
1090
      if ((double) p->blue < channel_statistics[BlueChannel].minima)
902
1091
        channel_statistics[BlueChannel].minima=(double) p->blue;
903
1092
      if ((double) p->blue > channel_statistics[BlueChannel].maxima)
904
1093
        channel_statistics[BlueChannel].maxima=(double) p->blue;
905
1094
      channel_statistics[BlueChannel].mean+=p->blue;
906
 
      channel_statistics[BlueChannel].standard_deviation+=(double)
 
1095
      channel_statistics[BlueChannel].standard_deviation+=(double) p->blue*
 
1096
        p->blue;
 
1097
      channel_statistics[BlueChannel].kurtosis+=(double) p->blue*p->blue*
907
1098
        p->blue*p->blue;
 
1099
      channel_statistics[BlueChannel].skewness+=(double) p->blue*p->blue*
 
1100
        p->blue;
908
1101
      if ((double) p->opacity < channel_statistics[OpacityChannel].minima)
909
1102
        channel_statistics[OpacityChannel].minima=(double) p->opacity;
910
1103
      if ((double) p->opacity > channel_statistics[OpacityChannel].maxima)
912
1105
      channel_statistics[OpacityChannel].mean+=p->opacity;
913
1106
      channel_statistics[OpacityChannel].standard_deviation+=(double)
914
1107
        p->opacity*p->opacity;
 
1108
      channel_statistics[OpacityChannel].kurtosis+=(double) p->opacity*
 
1109
        p->opacity*p->opacity*p->opacity;
 
1110
      channel_statistics[OpacityChannel].skewness+=(double) p->opacity*
 
1111
        p->opacity*p->opacity;
915
1112
      if (image->colorspace == CMYKColorspace)
916
1113
        {
917
1114
          if ((double) indexes[x] < channel_statistics[BlackChannel].minima)
921
1118
          channel_statistics[BlackChannel].mean+=indexes[x];
922
1119
          channel_statistics[BlackChannel].standard_deviation+=(double)
923
1120
            indexes[x]*indexes[x];
 
1121
          channel_statistics[BlackChannel].kurtosis+=(double) indexes[x]*
 
1122
            indexes[x]*indexes[x]*indexes[x];
 
1123
          channel_statistics[BlackChannel].skewness+=(double) indexes[x]*
 
1124
            indexes[x]*indexes[x];
924
1125
        }
925
1126
      x++;
926
1127
      p++;
931
1132
  {
932
1133
    channel_statistics[i].mean/=area;
933
1134
    channel_statistics[i].standard_deviation/=area;
 
1135
    channel_statistics[i].kurtosis/=area;
 
1136
    channel_statistics[i].skewness/=area;
934
1137
  }
935
1138
  for (i=0; i < AllChannels; i++)
936
1139
  {
937
 
    channel_statistics[AllChannels].depth=(unsigned long) MagickMax((double) 
 
1140
    channel_statistics[AllChannels].depth=(unsigned long) MagickMax((double)
938
1141
      channel_statistics[AllChannels].depth,(double)
939
1142
      channel_statistics[i].depth);
940
1143
    channel_statistics[AllChannels].minima=MagickMin(
944
1147
    channel_statistics[AllChannels].mean+=channel_statistics[i].mean;
945
1148
    channel_statistics[AllChannels].standard_deviation+=
946
1149
      channel_statistics[i].standard_deviation;
 
1150
    channel_statistics[AllChannels].kurtosis+=channel_statistics[i].kurtosis;
 
1151
    channel_statistics[AllChannels].skewness+=channel_statistics[i].skewness;
947
1152
  }
948
1153
  channels=4;
949
1154
  if (image->colorspace == CMYKColorspace)
950
1155
    channels++;
951
1156
  channel_statistics[AllChannels].mean/=channels;
952
1157
  channel_statistics[AllChannels].standard_deviation/=channels;
 
1158
  channel_statistics[AllChannels].kurtosis/=channels;
 
1159
  channel_statistics[AllChannels].skewness/=channels;
953
1160
  for (i=0; i <= AllChannels; i++)
 
1161
  {
 
1162
    sum_squares=0.0;
 
1163
    sum_squares=channel_statistics[i].standard_deviation;
 
1164
    sum_cubes=0.0;
 
1165
    sum_cubes=channel_statistics[i].skewness;
954
1166
    channel_statistics[i].standard_deviation=sqrt(
955
1167
      channel_statistics[i].standard_deviation-
956
1168
       (channel_statistics[i].mean*channel_statistics[i].mean));
 
1169
    if (channel_statistics[i].standard_deviation == 0.0)
 
1170
      {
 
1171
        channel_statistics[i].kurtosis=0.0;
 
1172
        channel_statistics[i].skewness=0.0;
 
1173
      }
 
1174
    else
 
1175
      {
 
1176
        channel_statistics[i].skewness=(channel_statistics[i].skewness-
 
1177
          3.0*channel_statistics[i].mean*sum_squares+
 
1178
          2.0*channel_statistics[i].mean*channel_statistics[i].mean*
 
1179
          channel_statistics[i].mean)/
 
1180
          (channel_statistics[i].standard_deviation*
 
1181
           channel_statistics[i].standard_deviation*
 
1182
           channel_statistics[i].standard_deviation);
 
1183
        channel_statistics[i].kurtosis=(channel_statistics[i].kurtosis-
 
1184
          4.0*channel_statistics[i].mean*sum_cubes+
 
1185
          6.0*channel_statistics[i].mean*channel_statistics[i].mean*sum_squares-
 
1186
          3.0*channel_statistics[i].mean*channel_statistics[i].mean*
 
1187
          1.0*channel_statistics[i].mean*channel_statistics[i].mean)/
 
1188
          (channel_statistics[i].standard_deviation*
 
1189
           channel_statistics[i].standard_deviation*
 
1190
           channel_statistics[i].standard_deviation*
 
1191
           channel_statistics[i].standard_deviation)-3.0;
 
1192
      }
 
1193
  }
957
1194
  return(channel_statistics);
958
1195
}
959
1196
 
1046
1283
MagickExport MagickBooleanType SetImageChannelDepth(Image *image,
1047
1284
  const ChannelType channel,const unsigned long depth)
1048
1285
{
 
1286
  ExceptionInfo
 
1287
    *exception;
 
1288
 
1049
1289
  long
1050
1290
    y;
1051
1291
 
1053
1293
    status;
1054
1294
 
1055
1295
  QuantumAny
1056
 
    scale;
 
1296
    range;
1057
1297
 
1058
1298
  ViewInfo
1059
 
    **image_view;
1060
 
  
 
1299
    *image_view;
 
1300
 
1061
1301
  assert(image != (Image *) NULL);
1062
1302
  if (image->debug != MagickFalse)
1063
1303
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
1072
1312
    Scale pixels to desired depth.
1073
1313
  */
1074
1314
  status=MagickTrue;
1075
 
  scale=GetQuantumScale(depth);
1076
 
  image_view=AcquireCacheViewThreadSet(image);
 
1315
  range=GetQuantumRange(depth);
 
1316
  exception=(&image->exception);
 
1317
  image_view=AcquireCacheView(image);
1077
1318
#if defined(MAGICKCORE_OPENMP_SUPPORT)
1078
 
  #pragma omp parallel for schedule(static,64) shared(status)
 
1319
  #pragma omp parallel for schedule(dynamic,4) shared(status)
1079
1320
#endif
1080
1321
  for (y=0; y < (long) image->rows; y++)
1081
1322
  {
1083
1324
      *indexes;
1084
1325
 
1085
1326
    register long
1086
 
      id,
1087
1327
      x;
1088
1328
 
1089
1329
    register PixelPacket
1091
1331
 
1092
1332
    if (status == MagickFalse)
1093
1333
      continue;
1094
 
    id=GetCacheViewThreadId();
1095
 
    q=GetCacheViewPixels(image_view[id],0,y,image->columns,1);
 
1334
    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,
 
1335
      exception);
1096
1336
    if (q == (PixelPacket *) NULL)
1097
1337
      {
1098
1338
        status=MagickFalse;
1099
1339
        continue;
1100
1340
      }
1101
 
    indexes=GetCacheViewIndexes(image_view[id]);
 
1341
    indexes=GetCacheViewAuthenticIndexQueue(image_view);
1102
1342
    for (x=0; x < (long) image->columns; x++)
1103
1343
    {
1104
1344
      if ((channel & RedChannel) != 0)
1105
 
        q->red=ScaleAnyToQuantum(ScaleQuantumToAny(q->red,depth,scale),depth,
1106
 
          scale);
 
1345
        q->red=ScaleAnyToQuantum(ScaleQuantumToAny(q->red,range),range);
1107
1346
      if ((channel & GreenChannel) != 0)
1108
 
        q->green=ScaleAnyToQuantum(ScaleQuantumToAny(q->green,depth,scale),
1109
 
          depth,scale);
 
1347
        q->green=ScaleAnyToQuantum(ScaleQuantumToAny(q->green,range),range);
1110
1348
      if ((channel & BlueChannel) != 0)
1111
 
        q->blue=ScaleAnyToQuantum(ScaleQuantumToAny(q->blue,depth,scale),
1112
 
          depth,scale);
 
1349
        q->blue=ScaleAnyToQuantum(ScaleQuantumToAny(q->blue,range),range);
1113
1350
      if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1114
 
        q->opacity=ScaleAnyToQuantum(ScaleQuantumToAny(q->opacity,depth,scale),
1115
 
          depth,scale);
 
1351
        q->opacity=ScaleAnyToQuantum(ScaleQuantumToAny(q->opacity,range),range);
1116
1352
      if (((channel & IndexChannel) != 0) &&
1117
1353
          (image->colorspace == CMYKColorspace))
1118
 
        indexes[x]=ScaleAnyToQuantum(ScaleQuantumToAny(indexes[x],depth,scale),
1119
 
          depth,scale);
 
1354
        indexes[x]=ScaleAnyToQuantum(ScaleQuantumToAny(indexes[x],range),range);
1120
1355
      q++;
1121
1356
    }
1122
 
    if (SyncCacheViewPixels(image_view[id]) == MagickFalse)
 
1357
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1123
1358
      {
1124
1359
        status=MagickFalse;
1125
1360
        continue;
1126
1361
      }
1127
1362
  }
1128
 
  image_view=DestroyCacheViewThreadSet(image_view);
 
1363
  image_view=DestroyCacheView(image_view);
1129
1364
  if (image->storage_class == PseudoClass)
1130
1365
    {
1131
1366
      QuantumAny
1132
 
        scale;
 
1367
        range;
1133
1368
 
1134
1369
      register long
1135
1370
        i;
1138
1373
        *p;
1139
1374
 
1140
1375
      p=image->colormap;
1141
 
      scale=GetQuantumScale(depth);
 
1376
      range=GetQuantumRange(depth);
1142
1377
#if defined(MAGICKCORE_OPENMP_SUPPORT)
1143
 
      #pragma omp parallel for schedule(static,64) shared(status)
 
1378
  #pragma omp parallel for schedule(dynamic,4) shared(status)
1144
1379
#endif
1145
1380
      for (i=0; i < (long) image->colors; i++)
1146
1381
      {
1147
1382
        if ((channel & RedChannel) != 0)
1148
 
          p->red=ScaleAnyToQuantum(ScaleQuantumToAny(p->red,depth,
1149
 
            scale),depth,scale);
 
1383
          p->red=ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),range);
1150
1384
        if ((channel & GreenChannel) != 0)
1151
 
          p->green=ScaleAnyToQuantum(ScaleQuantumToAny(p->green,depth,
1152
 
            scale),depth,scale);
 
1385
          p->green=ScaleAnyToQuantum(ScaleQuantumToAny(p->green,range),range);
1153
1386
        if ((channel & BlueChannel) != 0)
1154
 
          p->blue=ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,depth,
1155
 
            scale),depth,scale);
 
1387
          p->blue=ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,range),range);
1156
1388
        if ((channel & OpacityChannel) != 0)
1157
 
          p->opacity=ScaleAnyToQuantum(ScaleQuantumToAny(p->opacity,depth,
1158
 
            scale),depth,scale);
 
1389
          p->opacity=ScaleAnyToQuantum(ScaleQuantumToAny(p->opacity,range),
 
1390
            range);
1159
1391
        p++;
1160
1392
      }
1161
1393
    }