MagickCore  6.7.5
statistic.c
Go to the documentation of this file.
00001 /*
00002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00003 %                                                                             %
00004 %                                                                             %
00005 %                                                                             %
00006 %        SSSSS  TTTTT   AAA   TTTTT  IIIII  SSSSS  TTTTT  IIIII   CCCC        %
00007 %        SS       T    A   A    T      I    SS       T      I    C            %
00008 %         SSS     T    AAAAA    T      I     SSS     T      I    C            %
00009 %           SS    T    A   A    T      I       SS    T      I    C            %
00010 %        SSSSS    T    A   A    T    IIIII  SSSSS    T    IIIII   CCCC        %
00011 %                                                                             %
00012 %                                                                             %
00013 %                     MagickCore Image Statistical Methods                    %
00014 %                                                                             %
00015 %                              Software Design                                %
00016 %                                John Cristy                                  %
00017 %                                 July 1992                                   %
00018 %                                                                             %
00019 %                                                                             %
00020 %  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
00021 %  dedicated to making software imaging solutions freely available.           %
00022 %                                                                             %
00023 %  You may not use this file except in compliance with the License.  You may  %
00024 %  obtain a copy of the License at                                            %
00025 %                                                                             %
00026 %    http://www.imagemagick.org/script/license.php                            %
00027 %                                                                             %
00028 %  Unless required by applicable law or agreed to in writing, software        %
00029 %  distributed under the License is distributed on an "AS IS" BASIS,          %
00030 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
00031 %  See the License for the specific language governing permissions and        %
00032 %  limitations under the License.                                             %
00033 %                                                                             %
00034 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00035 %
00036 %
00037 %
00038 */
00039 
00040 /*
00041   Include declarations.
00042 */
00043 #include "MagickCore/studio.h"
00044 #include "MagickCore/property.h"
00045 #include "MagickCore/animate.h"
00046 #include "MagickCore/blob.h"
00047 #include "MagickCore/blob-private.h"
00048 #include "MagickCore/cache.h"
00049 #include "MagickCore/cache-private.h"
00050 #include "MagickCore/cache-view.h"
00051 #include "MagickCore/client.h"
00052 #include "MagickCore/color.h"
00053 #include "MagickCore/color-private.h"
00054 #include "MagickCore/colorspace.h"
00055 #include "MagickCore/colorspace-private.h"
00056 #include "MagickCore/composite.h"
00057 #include "MagickCore/composite-private.h"
00058 #include "MagickCore/compress.h"
00059 #include "MagickCore/constitute.h"
00060 #include "MagickCore/display.h"
00061 #include "MagickCore/draw.h"
00062 #include "MagickCore/enhance.h"
00063 #include "MagickCore/exception.h"
00064 #include "MagickCore/exception-private.h"
00065 #include "MagickCore/gem.h"
00066 #include "MagickCore/gem-private.h"
00067 #include "MagickCore/geometry.h"
00068 #include "MagickCore/list.h"
00069 #include "MagickCore/image-private.h"
00070 #include "MagickCore/magic.h"
00071 #include "MagickCore/magick.h"
00072 #include "MagickCore/memory_.h"
00073 #include "MagickCore/module.h"
00074 #include "MagickCore/monitor.h"
00075 #include "MagickCore/monitor-private.h"
00076 #include "MagickCore/option.h"
00077 #include "MagickCore/paint.h"
00078 #include "MagickCore/pixel-accessor.h"
00079 #include "MagickCore/profile.h"
00080 #include "MagickCore/quantize.h"
00081 #include "MagickCore/quantum-private.h"
00082 #include "MagickCore/random_.h"
00083 #include "MagickCore/random-private.h"
00084 #include "MagickCore/segment.h"
00085 #include "MagickCore/semaphore.h"
00086 #include "MagickCore/signature-private.h"
00087 #include "MagickCore/statistic.h"
00088 #include "MagickCore/string_.h"
00089 #include "MagickCore/thread-private.h"
00090 #include "MagickCore/timer.h"
00091 #include "MagickCore/utility.h"
00092 #include "MagickCore/version.h"
00093 
00094 /*
00095 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00096 %                                                                             %
00097 %                                                                             %
00098 %                                                                             %
00099 %     E v a l u a t e I m a g e                                               %
00100 %                                                                             %
00101 %                                                                             %
00102 %                                                                             %
00103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00104 %
00105 %  EvaluateImage() applies a value to the image with an arithmetic, relational,
00106 %  or logical operator to an image. Use these operations to lighten or darken
00107 %  an image, to increase or decrease contrast in an image, or to produce the
00108 %  "negative" of an image.
00109 %
00110 %  The format of the EvaluateImage method is:
00111 %
00112 %      MagickBooleanType EvaluateImage(Image *image,
00113 %        const MagickEvaluateOperator op,const double value,
00114 %        ExceptionInfo *exception)
00115 %      MagickBooleanType EvaluateImages(Image *images,
00116 %        const MagickEvaluateOperator op,const double value,
00117 %        ExceptionInfo *exception)
00118 %
00119 %  A description of each parameter follows:
00120 %
00121 %    o image: the image.
00122 %
00123 %    o op: A channel op.
00124 %
00125 %    o value: A value value.
00126 %
00127 %    o exception: return any errors or warnings in this structure.
00128 %
00129 */
00130 
00131 typedef struct _PixelChannels
00132 {
00133   MagickRealType
00134     channel[CompositePixelChannel];
00135 } PixelChannels;
00136 
00137 static PixelChannels **DestroyPixelThreadSet(PixelChannels **pixels)
00138 {
00139   register ssize_t
00140     i;
00141 
00142   assert(pixels != (PixelChannels **) NULL);
00143   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
00144     if (pixels[i] != (PixelChannels *) NULL)
00145       pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]);
00146   pixels=(PixelChannels **) RelinquishMagickMemory(pixels);
00147   return(pixels);
00148 }
00149 
00150 static PixelChannels **AcquirePixelThreadSet(const Image *image,
00151   const size_t number_images)
00152 {
00153   register ssize_t
00154     i;
00155 
00156   PixelChannels
00157     **pixels;
00158 
00159   size_t
00160     length,
00161     number_threads;
00162 
00163   number_threads=GetOpenMPMaximumThreads();
00164   pixels=(PixelChannels **) AcquireQuantumMemory(number_threads,
00165     sizeof(*pixels));
00166   if (pixels == (PixelChannels **) NULL)
00167     return((PixelChannels **) NULL);
00168   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
00169   for (i=0; i < (ssize_t) number_threads; i++)
00170   {
00171     register ssize_t
00172       j;
00173 
00174     length=image->columns;
00175     if (length < number_images)
00176       length=number_images;
00177     pixels[i]=(PixelChannels *) AcquireQuantumMemory(length,sizeof(**pixels));
00178     if (pixels[i] == (PixelChannels *) NULL)
00179       return(DestroyPixelThreadSet(pixels));
00180     for (j=0; j < (ssize_t) length; j++)
00181     {
00182       register ssize_t
00183         k;
00184 
00185       for (k=0; k < MaxPixelChannels; k++)
00186         pixels[i][j].channel[k]=0.0;
00187     }
00188   }
00189   return(pixels);
00190 }
00191 
00192 static inline double EvaluateMax(const double x,const double y)
00193 {
00194   if (x > y)
00195     return(x);
00196   return(y);
00197 }
00198 
00199 #if defined(__cplusplus) || defined(c_plusplus)
00200 extern "C" {
00201 #endif
00202 
00203 static int IntensityCompare(const void *x,const void *y)
00204 {
00205   const PixelChannels
00206     *color_1,
00207     *color_2;
00208 
00209   MagickRealType
00210     distance;
00211 
00212   register ssize_t
00213     i;
00214 
00215   color_1=(const PixelChannels *) x;
00216   color_2=(const PixelChannels *) y;
00217   distance=0.0;
00218   for (i=0; i < MaxPixelChannels; i++)
00219     distance+=color_1->channel[i]-(MagickRealType) color_2->channel[i];
00220   return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
00221 }
00222 
00223 #if defined(__cplusplus) || defined(c_plusplus)
00224 }
00225 #endif
00226 
00227 static inline double MagickMin(const double x,const double y)
00228 {
00229   if (x < y)
00230     return(x);
00231   return(y);
00232 }
00233 
00234 static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
00235   Quantum pixel,const MagickEvaluateOperator op,const MagickRealType value)
00236 {
00237   MagickRealType
00238     result;
00239 
00240   result=0.0;
00241   switch (op)
00242   {
00243     case UndefinedEvaluateOperator:
00244       break;
00245     case AbsEvaluateOperator:
00246     {
00247       result=(MagickRealType) fabs((double) (pixel+value));
00248       break;
00249     }
00250     case AddEvaluateOperator:
00251     {
00252       result=(MagickRealType) (pixel+value);
00253       break;
00254     }
00255     case AddModulusEvaluateOperator:
00256     {
00257       /*
00258         This returns a 'floored modulus' of the addition which is a positive
00259         result.  It differs from % or fmod() that returns a 'truncated modulus'
00260         result, where floor() is replaced by trunc() and could return a
00261         negative result (which is clipped).
00262       */
00263       result=pixel+value;
00264       result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
00265       break;
00266     }
00267     case AndEvaluateOperator:
00268     {
00269       result=(MagickRealType) ((size_t) pixel & (size_t) (value+0.5));
00270       break;
00271     }
00272     case CosineEvaluateOperator:
00273     {
00274       result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
00275         QuantumScale*pixel*value))+0.5));
00276       break;
00277     }
00278     case DivideEvaluateOperator:
00279     {
00280       result=pixel/(value == 0.0 ? 1.0 : value);
00281       break;
00282     }
00283     case ExponentialEvaluateOperator:
00284     {
00285       result=(MagickRealType) (QuantumRange*exp((double) (value*QuantumScale*
00286         pixel)));
00287       break;
00288     }
00289     case GaussianNoiseEvaluateOperator:
00290     {
00291       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
00292         GaussianNoise,value);
00293       break;
00294     }
00295     case ImpulseNoiseEvaluateOperator:
00296     {
00297       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
00298         ImpulseNoise,value);
00299       break;
00300     }
00301     case LaplacianNoiseEvaluateOperator:
00302     {
00303       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
00304         LaplacianNoise,value);
00305       break;
00306     }
00307     case LeftShiftEvaluateOperator:
00308     {
00309       result=(MagickRealType) ((size_t) pixel << (size_t) (value+0.5));
00310       break;
00311     }
00312     case LogEvaluateOperator:
00313     {
00314       result=(MagickRealType) (QuantumRange*log((double) (QuantumScale*value*
00315         pixel+1.0))/log((double) (value+1.0)));
00316       break;
00317     }
00318     case MaxEvaluateOperator:
00319     {
00320       result=(MagickRealType) EvaluateMax((double) pixel,value);
00321       break;
00322     }
00323     case MeanEvaluateOperator:
00324     {
00325       result=(MagickRealType) (pixel+value);
00326       break;
00327     }
00328     case MedianEvaluateOperator:
00329     {
00330       result=(MagickRealType) (pixel+value);
00331       break;
00332     }
00333     case MinEvaluateOperator:
00334     {
00335       result=(MagickRealType) MagickMin((double) pixel,value);
00336       break;
00337     }
00338     case MultiplicativeNoiseEvaluateOperator:
00339     {
00340       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
00341         MultiplicativeGaussianNoise,value);
00342       break;
00343     }
00344     case MultiplyEvaluateOperator:
00345     {
00346       result=(MagickRealType) (value*pixel);
00347       break;
00348     }
00349     case OrEvaluateOperator:
00350     {
00351       result=(MagickRealType) ((size_t) pixel | (size_t) (value+0.5));
00352       break;
00353     }
00354     case PoissonNoiseEvaluateOperator:
00355     {
00356       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
00357         PoissonNoise,value);
00358       break;
00359     }
00360     case PowEvaluateOperator:
00361     {
00362       result=(MagickRealType) (QuantumRange*pow((double) (QuantumScale*pixel),
00363         (double) value));
00364       break;
00365     }
00366     case RightShiftEvaluateOperator:
00367     {
00368       result=(MagickRealType) ((size_t) pixel >> (size_t) (value+0.5));
00369       break;
00370     }
00371     case SetEvaluateOperator:
00372     {
00373       result=value;
00374       break;
00375     }
00376     case SineEvaluateOperator:
00377     {
00378       result=(MagickRealType) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
00379         QuantumScale*pixel*value))+0.5));
00380       break;
00381     }
00382     case SubtractEvaluateOperator:
00383     {
00384       result=(MagickRealType) (pixel-value);
00385       break;
00386     }
00387     case SumEvaluateOperator:
00388     {
00389       result=(MagickRealType) (pixel+value);
00390       break;
00391     }
00392     case ThresholdEvaluateOperator:
00393     {
00394       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 :
00395         QuantumRange);
00396       break;
00397     }
00398     case ThresholdBlackEvaluateOperator:
00399     {
00400       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : pixel);
00401       break;
00402     }
00403     case ThresholdWhiteEvaluateOperator:
00404     {
00405       result=(MagickRealType) (((MagickRealType) pixel > value) ? QuantumRange :
00406         pixel);
00407       break;
00408     }
00409     case UniformNoiseEvaluateOperator:
00410     {
00411       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
00412         UniformNoise,value);
00413       break;
00414     }
00415     case XorEvaluateOperator:
00416     {
00417       result=(MagickRealType) ((size_t) pixel ^ (size_t) (value+0.5));
00418       break;
00419     }
00420   }
00421   return(result);
00422 }
00423 
00424 MagickExport Image *EvaluateImages(const Image *images,
00425   const MagickEvaluateOperator op,ExceptionInfo *exception)
00426 {
00427 #define EvaluateImageTag  "Evaluate/Image"
00428 
00429   CacheView
00430     *evaluate_view;
00431 
00432   const Image
00433     *next;
00434 
00435   Image
00436     *evaluate_image;
00437 
00438   MagickBooleanType
00439     status;
00440 
00441   MagickOffsetType
00442     progress;
00443 
00444   PixelChannels
00445     **restrict evaluate_pixels;
00446 
00447   RandomInfo
00448     **restrict random_info;
00449 
00450   size_t
00451     number_images;
00452 
00453   ssize_t
00454     y;
00455 
00456   /*
00457     Ensure the image are the same size.
00458   */
00459   assert(images != (Image *) NULL);
00460   assert(images->signature == MagickSignature);
00461   if (images->debug != MagickFalse)
00462     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
00463   assert(exception != (ExceptionInfo *) NULL);
00464   assert(exception->signature == MagickSignature);
00465   for (next=images; next != (Image *) NULL; next=GetNextImageInList(next))
00466     if ((next->columns != images->columns) || (next->rows != images->rows))
00467       {
00468         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00469           "ImageWidthsOrHeightsDiffer","`%s'",images->filename);
00470         return((Image *) NULL);
00471       }
00472   /*
00473     Initialize evaluate next attributes.
00474   */
00475   evaluate_image=CloneImage(images,images->columns,images->rows,MagickTrue,
00476     exception);
00477   if (evaluate_image == (Image *) NULL)
00478     return((Image *) NULL);
00479   if (SetImageStorageClass(evaluate_image,DirectClass,exception) == MagickFalse)
00480     {
00481       evaluate_image=DestroyImage(evaluate_image);
00482       return((Image *) NULL);
00483     }
00484   number_images=GetImageListLength(images);
00485   evaluate_pixels=AcquirePixelThreadSet(images,number_images);
00486   if (evaluate_pixels == (PixelChannels **) NULL)
00487     {
00488       evaluate_image=DestroyImage(evaluate_image);
00489       (void) ThrowMagickException(exception,GetMagickModule(),
00490         ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
00491       return((Image *) NULL);
00492     }
00493   /*
00494     Evaluate image pixels.
00495   */
00496   status=MagickTrue;
00497   progress=0;
00498   random_info=AcquireRandomInfoThreadSet();
00499   evaluate_view=AcquireCacheView(evaluate_image);
00500   if (op == MedianEvaluateOperator)
00501     {
00502 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00503       #pragma omp parallel for schedule(static) shared(progress,status)
00504 #endif
00505       for (y=0; y < (ssize_t) evaluate_image->rows; y++)
00506       {
00507         CacheView
00508           *image_view;
00509 
00510         const Image
00511           *next;
00512 
00513         const int
00514           id = GetOpenMPThreadId();
00515 
00516         register PixelChannels
00517           *evaluate_pixel;
00518 
00519         register Quantum
00520           *restrict q;
00521 
00522         register ssize_t
00523           x;
00524 
00525         if (status == MagickFalse)
00526           continue;
00527         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
00528           evaluate_image->columns,1,exception);
00529         if (q == (Quantum *) NULL)
00530           {
00531             status=MagickFalse;
00532             continue;
00533           }
00534         evaluate_pixel=evaluate_pixels[id];
00535         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
00536         {
00537           register ssize_t
00538             j,
00539             k;
00540 
00541           for (j=0; j < (ssize_t) number_images; j++)
00542             for (k=0; k < MaxPixelChannels; k++)
00543               evaluate_pixel[j].channel[k]=0.0;
00544           next=images;
00545           for (j=0; j < (ssize_t) number_images; j++)
00546           {
00547             register const Quantum
00548               *p;
00549 
00550             register ssize_t
00551               i;
00552 
00553             image_view=AcquireCacheView(next);
00554             p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
00555             if (p == (const Quantum *) NULL)
00556               {
00557                 image_view=DestroyCacheView(image_view);
00558                 break;
00559               }
00560             for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
00561             {
00562               PixelChannel
00563                 channel;
00564 
00565               PixelTrait
00566                 evaluate_traits,
00567                 traits;
00568 
00569               channel=GetPixelChannelMapChannel(evaluate_image,i);
00570               evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
00571               traits=GetPixelChannelMapTraits(next,channel);
00572               if ((traits == UndefinedPixelTrait) ||
00573                   (evaluate_traits == UndefinedPixelTrait))
00574                 continue;
00575               if ((evaluate_traits & UpdatePixelTrait) == 0)
00576                 continue;
00577               evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
00578                 random_info[id],GetPixelChannel(evaluate_image,channel,p),op,
00579                 evaluate_pixel[j].channel[i]);
00580             }
00581             image_view=DestroyCacheView(image_view);
00582             next=GetNextImageInList(next);
00583           }
00584           qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
00585             IntensityCompare);
00586           for (k=0; k < (ssize_t) GetPixelChannels(evaluate_image); k++)
00587             q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
00588           q+=GetPixelChannels(evaluate_image);
00589         }
00590         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
00591           status=MagickFalse;
00592         if (images->progress_monitor != (MagickProgressMonitor) NULL)
00593           {
00594             MagickBooleanType
00595               proceed;
00596 
00597 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
00598             #pragma omp critical (MagickCore_EvaluateImages)
00599 #endif
00600             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
00601               evaluate_image->rows);
00602             if (proceed == MagickFalse)
00603               status=MagickFalse;
00604           }
00605       }
00606     }
00607   else
00608     {
00609 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00610       #pragma omp parallel for schedule(static) shared(progress,status)
00611 #endif
00612       for (y=0; y < (ssize_t) evaluate_image->rows; y++)
00613       {
00614         CacheView
00615           *image_view;
00616 
00617         const Image
00618           *next;
00619 
00620         const int
00621           id = GetOpenMPThreadId();
00622 
00623         register ssize_t
00624           i,
00625           x;
00626 
00627         register PixelChannels
00628           *evaluate_pixel;
00629 
00630         register Quantum
00631           *restrict q;
00632 
00633         ssize_t
00634           j;
00635 
00636         if (status == MagickFalse)
00637           continue;
00638         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
00639           evaluate_image->columns,1,exception);
00640         if (q == (Quantum *) NULL)
00641           {
00642             status=MagickFalse;
00643             continue;
00644           }
00645         evaluate_pixel=evaluate_pixels[id];
00646         for (j=0; j < (ssize_t) evaluate_image->columns; j++)
00647           for (i=0; i < MaxPixelChannels; i++)
00648             evaluate_pixel[j].channel[i]=0.0;
00649         next=images;
00650         for (j=0; j < (ssize_t) number_images; j++)
00651         {
00652           register const Quantum
00653             *p;
00654 
00655           image_view=AcquireCacheView(next);
00656           p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
00657           if (p == (const Quantum *) NULL)
00658             {
00659               image_view=DestroyCacheView(image_view);
00660               break;
00661             }
00662           for (x=0; x < (ssize_t) next->columns; x++)
00663           {
00664             register ssize_t
00665               i;
00666 
00667             if (GetPixelMask(next,p) != 0)
00668               {
00669                 p+=GetPixelChannels(next);
00670                 continue;
00671               }
00672             for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
00673             {
00674               PixelChannel
00675                 channel;
00676 
00677               PixelTrait
00678                 evaluate_traits,
00679                 traits;
00680 
00681               channel=GetPixelChannelMapChannel(evaluate_image,i);
00682               traits=GetPixelChannelMapTraits(next,channel);
00683               evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
00684               if ((traits == UndefinedPixelTrait) ||
00685                   (evaluate_traits == UndefinedPixelTrait))
00686                 continue;
00687               if ((traits & UpdatePixelTrait) == 0)
00688                 continue;
00689               evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
00690                 random_info[id],GetPixelChannel(evaluate_image,channel,p),j ==
00691                 0 ? AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
00692             }
00693             p+=GetPixelChannels(next);
00694           }
00695           image_view=DestroyCacheView(image_view);
00696           next=GetNextImageInList(next);
00697         }
00698         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
00699         {
00700           register ssize_t
00701              i;
00702 
00703           switch (op)
00704           {
00705             case MeanEvaluateOperator:
00706             {
00707               for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
00708                 evaluate_pixel[x].channel[i]/=(MagickRealType) number_images;
00709               break;
00710             }
00711             case MultiplyEvaluateOperator:
00712             {
00713               for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
00714               {
00715                 register ssize_t
00716                   j;
00717 
00718                 for (j=0; j < (ssize_t) (number_images-1); j++)
00719                   evaluate_pixel[x].channel[i]*=QuantumScale;
00720               }
00721               break;
00722             }
00723             default:
00724               break;
00725           }
00726         }
00727         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
00728         {
00729           register ssize_t
00730             i;
00731 
00732           if (GetPixelMask(evaluate_image,q) != 0)
00733             {
00734               q+=GetPixelChannels(evaluate_image);
00735               continue;
00736             }
00737           for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
00738           {
00739             PixelChannel
00740               channel;
00741 
00742             PixelTrait
00743               traits;
00744 
00745             channel=GetPixelChannelMapChannel(evaluate_image,i);
00746             traits=GetPixelChannelMapTraits(evaluate_image,channel);
00747             if (traits == UndefinedPixelTrait)
00748               continue;
00749             if ((traits & UpdatePixelTrait) == 0)
00750               continue;
00751             q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
00752           }
00753           q+=GetPixelChannels(evaluate_image);
00754         }
00755         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
00756           status=MagickFalse;
00757         if (images->progress_monitor != (MagickProgressMonitor) NULL)
00758           {
00759             MagickBooleanType
00760               proceed;
00761 
00762 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
00763             #pragma omp critical (MagickCore_EvaluateImages)
00764 #endif
00765             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
00766               evaluate_image->rows);
00767             if (proceed == MagickFalse)
00768               status=MagickFalse;
00769           }
00770       }
00771     }
00772   evaluate_view=DestroyCacheView(evaluate_view);
00773   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
00774   random_info=DestroyRandomInfoThreadSet(random_info);
00775   if (status == MagickFalse)
00776     evaluate_image=DestroyImage(evaluate_image);
00777   return(evaluate_image);
00778 }
00779 
00780 MagickExport MagickBooleanType EvaluateImage(Image *image,
00781   const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
00782 {
00783   CacheView
00784     *image_view;
00785 
00786   MagickBooleanType
00787     status;
00788 
00789   MagickOffsetType
00790     progress;
00791 
00792   RandomInfo
00793     **restrict random_info;
00794 
00795   ssize_t
00796     y;
00797 
00798   assert(image != (Image *) NULL);
00799   assert(image->signature == MagickSignature);
00800   if (image->debug != MagickFalse)
00801     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00802   assert(exception != (ExceptionInfo *) NULL);
00803   assert(exception->signature == MagickSignature);
00804   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
00805     return(MagickFalse);
00806   status=MagickTrue;
00807   progress=0;
00808   random_info=AcquireRandomInfoThreadSet();
00809   image_view=AcquireCacheView(image);
00810 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00811   #pragma omp parallel for schedule(static,4) shared(progress,status)
00812 #endif
00813   for (y=0; y < (ssize_t) image->rows; y++)
00814   {
00815     const int
00816       id = GetOpenMPThreadId();
00817 
00818     register Quantum
00819       *restrict q;
00820 
00821     register ssize_t
00822       x;
00823 
00824     if (status == MagickFalse)
00825       continue;
00826     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
00827     if (q == (Quantum *) NULL)
00828       {
00829         status=MagickFalse;
00830         continue;
00831       }
00832     for (x=0; x < (ssize_t) image->columns; x++)
00833     {
00834       register ssize_t
00835         i;
00836 
00837       if (GetPixelMask(image,q) != 0)
00838         {
00839           q+=GetPixelChannels(image);
00840           continue;
00841         }
00842       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
00843       {
00844         PixelChannel
00845           channel;
00846 
00847         PixelTrait
00848           traits;
00849 
00850         channel=GetPixelChannelMapChannel(image,i);
00851         traits=GetPixelChannelMapTraits(image,channel);
00852         if (traits == UndefinedPixelTrait)
00853           continue;
00854         q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op,
00855           value));
00856       }
00857       q+=GetPixelChannels(image);
00858     }
00859     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
00860       status=MagickFalse;
00861     if (image->progress_monitor != (MagickProgressMonitor) NULL)
00862       {
00863         MagickBooleanType
00864           proceed;
00865 
00866 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00867   #pragma omp critical (MagickCore_EvaluateImage)
00868 #endif
00869         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
00870         if (proceed == MagickFalse)
00871           status=MagickFalse;
00872       }
00873   }
00874   image_view=DestroyCacheView(image_view);
00875   random_info=DestroyRandomInfoThreadSet(random_info);
00876   return(status);
00877 }
00878 
00879 /*
00880 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00881 %                                                                             %
00882 %                                                                             %
00883 %                                                                             %
00884 %     F u n c t i o n I m a g e                                               %
00885 %                                                                             %
00886 %                                                                             %
00887 %                                                                             %
00888 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00889 %
00890 %  FunctionImage() applies a value to the image with an arithmetic, relational,
00891 %  or logical operator to an image. Use these operations to lighten or darken
00892 %  an image, to increase or decrease contrast in an image, or to produce the
00893 %  "negative" of an image.
00894 %
00895 %  The format of the FunctionImage method is:
00896 %
00897 %      MagickBooleanType FunctionImage(Image *image,
00898 %        const MagickFunction function,const ssize_t number_parameters,
00899 %        const double *parameters,ExceptionInfo *exception)
00900 %
00901 %  A description of each parameter follows:
00902 %
00903 %    o image: the image.
00904 %
00905 %    o function: A channel function.
00906 %
00907 %    o parameters: one or more parameters.
00908 %
00909 %    o exception: return any errors or warnings in this structure.
00910 %
00911 */
00912 
00913 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
00914   const size_t number_parameters,const double *parameters,
00915   ExceptionInfo *exception)
00916 {
00917   MagickRealType
00918     result;
00919 
00920   register ssize_t
00921     i;
00922 
00923   (void) exception;
00924   result=0.0;
00925   switch (function)
00926   {
00927     case PolynomialFunction:
00928     {
00929       /*
00930         Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
00931         c1*x^2+c2*x+c3).
00932       */
00933       result=0.0;
00934       for (i=0; i < (ssize_t) number_parameters; i++)
00935         result=result*QuantumScale*pixel+parameters[i];
00936       result*=QuantumRange;
00937       break;
00938     }
00939     case SinusoidFunction:
00940     {
00941       MagickRealType
00942         amplitude,
00943         bias,
00944         frequency,
00945         phase;
00946 
00947       /*
00948         Sinusoid: frequency, phase, amplitude, bias.
00949       */
00950       frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
00951       phase=(number_parameters >= 2) ? parameters[1] : 0.0;
00952       amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
00953       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
00954       result=(MagickRealType) (QuantumRange*(amplitude*sin((double) (2.0*
00955         MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
00956       break;
00957     }
00958     case ArcsinFunction:
00959     {
00960       MagickRealType
00961         bias,
00962         center,
00963         range,
00964         width;
00965 
00966       /*
00967         Arcsin (peged at range limits for invalid results): width, center,
00968         range, and bias.
00969       */
00970       width=(number_parameters >= 1) ? parameters[0] : 1.0;
00971       center=(number_parameters >= 2) ? parameters[1] : 0.5;
00972       range=(number_parameters >= 3) ? parameters[2] : 1.0;
00973       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
00974       result=2.0/width*(QuantumScale*pixel-center);
00975       if ( result <= -1.0 )
00976         result=bias-range/2.0;
00977       else
00978         if (result >= 1.0)
00979           result=bias+range/2.0;
00980         else
00981           result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
00982       result*=QuantumRange;
00983       break;
00984     }
00985     case ArctanFunction:
00986     {
00987       MagickRealType
00988         center,
00989         bias,
00990         range,
00991         slope;
00992 
00993       /*
00994         Arctan: slope, center, range, and bias.
00995       */
00996       slope=(number_parameters >= 1) ? parameters[0] : 1.0;
00997       center=(number_parameters >= 2) ? parameters[1] : 0.5;
00998       range=(number_parameters >= 3) ? parameters[2] : 1.0;
00999       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
01000       result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
01001       result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
01002         result)+bias));
01003       break;
01004     }
01005     case UndefinedFunction:
01006       break;
01007   }
01008   return(ClampToQuantum(result));
01009 }
01010 
01011 MagickExport MagickBooleanType FunctionImage(Image *image,
01012   const MagickFunction function,const size_t number_parameters,
01013   const double *parameters,ExceptionInfo *exception)
01014 {
01015 #define FunctionImageTag  "Function/Image "
01016 
01017   CacheView
01018     *image_view;
01019 
01020   MagickBooleanType
01021     status;
01022 
01023   MagickOffsetType
01024     progress;
01025 
01026   ssize_t
01027     y;
01028 
01029   assert(image != (Image *) NULL);
01030   assert(image->signature == MagickSignature);
01031   if (image->debug != MagickFalse)
01032     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01033   assert(exception != (ExceptionInfo *) NULL);
01034   assert(exception->signature == MagickSignature);
01035   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
01036     return(MagickFalse);
01037   status=MagickTrue;
01038   progress=0;
01039   image_view=AcquireCacheView(image);
01040 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01041   #pragma omp parallel for schedule(static,4) shared(progress,status)
01042 #endif
01043   for (y=0; y < (ssize_t) image->rows; y++)
01044   {
01045     register Quantum
01046       *restrict q;
01047 
01048     register ssize_t
01049       x;
01050 
01051     if (status == MagickFalse)
01052       continue;
01053     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
01054     if (q == (Quantum *) NULL)
01055       {
01056         status=MagickFalse;
01057         continue;
01058       }
01059     for (x=0; x < (ssize_t) image->columns; x++)
01060     {
01061       register ssize_t
01062         i;
01063 
01064       if (GetPixelMask(image,q) != 0)
01065         {
01066           q+=GetPixelChannels(image);
01067           continue;
01068         }
01069       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
01070       {
01071         PixelChannel
01072           channel;
01073 
01074         PixelTrait
01075           traits;
01076 
01077         channel=GetPixelChannelMapChannel(image,i);
01078         traits=GetPixelChannelMapTraits(image,channel);
01079         if (traits == UndefinedPixelTrait)
01080           continue;
01081         if ((traits & UpdatePixelTrait) == 0)
01082           continue;
01083         q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
01084           exception);
01085       }
01086       q+=GetPixelChannels(image);
01087     }
01088     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
01089       status=MagickFalse;
01090     if (image->progress_monitor != (MagickProgressMonitor) NULL)
01091       {
01092         MagickBooleanType
01093           proceed;
01094 
01095 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01096   #pragma omp critical (MagickCore_FunctionImage)
01097 #endif
01098         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
01099         if (proceed == MagickFalse)
01100           status=MagickFalse;
01101       }
01102   }
01103   image_view=DestroyCacheView(image_view);
01104   return(status);
01105 }
01106 
01107 /*
01108 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01109 %                                                                             %
01110 %                                                                             %
01111 %                                                                             %
01112 %   G e t I m a g e E x t r e m a                                             %
01113 %                                                                             %
01114 %                                                                             %
01115 %                                                                             %
01116 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01117 %
01118 %  GetImageExtrema() returns the extrema of one or more image channels.
01119 %
01120 %  The format of the GetImageExtrema method is:
01121 %
01122 %      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
01123 %        size_t *maxima,ExceptionInfo *exception)
01124 %
01125 %  A description of each parameter follows:
01126 %
01127 %    o image: the image.
01128 %
01129 %    o minima: the minimum value in the channel.
01130 %
01131 %    o maxima: the maximum value in the channel.
01132 %
01133 %    o exception: return any errors or warnings in this structure.
01134 %
01135 */
01136 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
01137   size_t *minima,size_t *maxima,ExceptionInfo *exception)
01138 {
01139   double
01140     max,
01141     min;
01142 
01143   MagickBooleanType
01144     status;
01145 
01146   assert(image != (Image *) NULL);
01147   assert(image->signature == MagickSignature);
01148   if (image->debug != MagickFalse)
01149     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01150   status=GetImageRange(image,&min,&max,exception);
01151   *minima=(size_t) ceil(min-0.5);
01152   *maxima=(size_t) floor(max+0.5);
01153   return(status);
01154 }
01155 
01156 /*
01157 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01158 %                                                                             %
01159 %                                                                             %
01160 %                                                                             %
01161 %   G e t I m a g e M e a n                                                   %
01162 %                                                                             %
01163 %                                                                             %
01164 %                                                                             %
01165 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01166 %
01167 %  GetImageMean() returns the mean and standard deviation of one or more
01168 %  image channels.
01169 %
01170 %  The format of the GetImageMean method is:
01171 %
01172 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
01173 %        double *standard_deviation,ExceptionInfo *exception)
01174 %
01175 %  A description of each parameter follows:
01176 %
01177 %    o image: the image.
01178 %
01179 %    o mean: the average value in the channel.
01180 %
01181 %    o standard_deviation: the standard deviation of the channel.
01182 %
01183 %    o exception: return any errors or warnings in this structure.
01184 %
01185 */
01186 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
01187   double *standard_deviation,ExceptionInfo *exception)
01188 {
01189   ChannelStatistics
01190     *channel_statistics;
01191 
01192   register ssize_t
01193     i;
01194 
01195   size_t
01196     area;
01197 
01198   assert(image != (Image *) NULL);
01199   assert(image->signature == MagickSignature);
01200   if (image->debug != MagickFalse)
01201     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01202   channel_statistics=GetImageStatistics(image,exception);
01203   if (channel_statistics == (ChannelStatistics *) NULL)
01204     return(MagickFalse);
01205   area=0;
01206   channel_statistics[CompositePixelChannel].mean=0.0;
01207   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
01208   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
01209   {
01210     PixelChannel
01211       channel;
01212 
01213     PixelTrait
01214       traits;
01215 
01216     channel=GetPixelChannelMapChannel(image,i);
01217     traits=GetPixelChannelMapTraits(image,channel);
01218     if (traits == UndefinedPixelTrait)
01219       continue;
01220     if ((traits & UpdatePixelTrait) == 0)
01221       continue;
01222     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
01223     channel_statistics[CompositePixelChannel].standard_deviation+=
01224       channel_statistics[i].variance-channel_statistics[i].mean*
01225       channel_statistics[i].mean;
01226     area++;
01227   }
01228   channel_statistics[CompositePixelChannel].mean/=area;
01229   channel_statistics[CompositePixelChannel].standard_deviation=
01230     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
01231   *mean=channel_statistics[CompositePixelChannel].mean;
01232   *standard_deviation=
01233     channel_statistics[CompositePixelChannel].standard_deviation;
01234   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
01235     channel_statistics);
01236   return(MagickTrue);
01237 }
01238 
01239 /*
01240 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01241 %                                                                             %
01242 %                                                                             %
01243 %                                                                             %
01244 %   G e t I m a g e K u r t o s i s                                           %
01245 %                                                                             %
01246 %                                                                             %
01247 %                                                                             %
01248 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01249 %
01250 %  GetImageKurtosis() returns the kurtosis and skewness of one or more
01251 %  image channels.
01252 %
01253 %  The format of the GetImageKurtosis method is:
01254 %
01255 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
01256 %        double *skewness,ExceptionInfo *exception)
01257 %
01258 %  A description of each parameter follows:
01259 %
01260 %    o image: the image.
01261 %
01262 %    o kurtosis: the kurtosis of the channel.
01263 %
01264 %    o skewness: the skewness of the channel.
01265 %
01266 %    o exception: return any errors or warnings in this structure.
01267 %
01268 */
01269 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
01270   double *kurtosis,double *skewness,ExceptionInfo *exception)
01271 {
01272   CacheView
01273     *image_view;
01274 
01275   double
01276     area,
01277     mean,
01278     standard_deviation,
01279     sum_squares,
01280     sum_cubes,
01281     sum_fourth_power;
01282 
01283   MagickBooleanType
01284     status;
01285 
01286   ssize_t
01287     y;
01288 
01289   assert(image != (Image *) NULL);
01290   assert(image->signature == MagickSignature);
01291   if (image->debug != MagickFalse)
01292     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01293   status=MagickTrue;
01294   *kurtosis=0.0;
01295   *skewness=0.0;
01296   area=0.0;
01297   mean=0.0;
01298   standard_deviation=0.0;
01299   sum_squares=0.0;
01300   sum_cubes=0.0;
01301   sum_fourth_power=0.0;
01302   image_view=AcquireCacheView(image);
01303 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01304   #pragma omp parallel for schedule(static) shared(status)
01305 #endif
01306   for (y=0; y < (ssize_t) image->rows; y++)
01307   {
01308     register const Quantum
01309       *restrict p;
01310 
01311     register ssize_t
01312       x;
01313 
01314     if (status == MagickFalse)
01315       continue;
01316     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
01317     if (p == (const Quantum *) NULL)
01318       {
01319         status=MagickFalse;
01320         continue;
01321       }
01322     for (x=0; x < (ssize_t) image->columns; x++)
01323     {
01324       register ssize_t
01325         i;
01326 
01327       if (GetPixelMask(image,p) != 0)
01328         {
01329           p+=GetPixelChannels(image);
01330           continue;
01331         }
01332       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
01333       {
01334         PixelChannel
01335           channel;
01336 
01337         PixelTrait
01338           traits;
01339 
01340         channel=GetPixelChannelMapChannel(image,i);
01341         traits=GetPixelChannelMapTraits(image,channel);
01342         if (traits == UndefinedPixelTrait)
01343           continue;
01344         if ((traits & UpdatePixelTrait) == 0)
01345           continue;
01346 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01347         #pragma omp critical (MagickCore_GetImageKurtosis)
01348 #endif
01349         {
01350           mean+=p[i];
01351           sum_squares+=(double) p[i]*p[i];
01352           sum_cubes+=(double) p[i]*p[i]*p[i];
01353           sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
01354           area++;
01355         }
01356       }
01357       p+=GetPixelChannels(image);
01358     }
01359   }
01360   image_view=DestroyCacheView(image_view);
01361   if (area != 0.0)
01362     {
01363       mean/=area;
01364       sum_squares/=area;
01365       sum_cubes/=area;
01366       sum_fourth_power/=area;
01367     }
01368   standard_deviation=sqrt(sum_squares-(mean*mean));
01369   if (standard_deviation != 0.0)
01370     {
01371       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
01372         3.0*mean*mean*mean*mean;
01373       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
01374         standard_deviation;
01375       *kurtosis-=3.0;
01376       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
01377       *skewness/=standard_deviation*standard_deviation*standard_deviation;
01378     }
01379   return(status);
01380 }
01381 
01382 /*
01383 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01384 %                                                                             %
01385 %                                                                             %
01386 %                                                                             %
01387 %   G e t I m a g e R a n g e                                                 %
01388 %                                                                             %
01389 %                                                                             %
01390 %                                                                             %
01391 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01392 %
01393 %  GetImageRange() returns the range of one or more image channels.
01394 %
01395 %  The format of the GetImageRange method is:
01396 %
01397 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
01398 %        double *maxima,ExceptionInfo *exception)
01399 %
01400 %  A description of each parameter follows:
01401 %
01402 %    o image: the image.
01403 %
01404 %    o minima: the minimum value in the channel.
01405 %
01406 %    o maxima: the maximum value in the channel.
01407 %
01408 %    o exception: return any errors or warnings in this structure.
01409 %
01410 */
01411 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
01412   double *maxima,ExceptionInfo *exception)
01413 {
01414   CacheView
01415     *image_view;
01416 
01417   MagickBooleanType
01418     status;
01419 
01420   ssize_t
01421     y;
01422 
01423   assert(image != (Image *) NULL);
01424   assert(image->signature == MagickSignature);
01425   if (image->debug != MagickFalse)
01426     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01427   status=MagickTrue;
01428   *maxima=(-MagickHuge);
01429   *minima=MagickHuge;
01430   image_view=AcquireCacheView(image);
01431 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01432   #pragma omp parallel for schedule(static) shared(status)
01433 #endif
01434   for (y=0; y < (ssize_t) image->rows; y++)
01435   {
01436     register const Quantum
01437       *restrict p;
01438 
01439     register ssize_t
01440       x;
01441 
01442     if (status == MagickFalse)
01443       continue;
01444     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
01445     if (p == (const Quantum *) NULL)
01446       {
01447         status=MagickFalse;
01448         continue;
01449       }
01450     for (x=0; x < (ssize_t) image->columns; x++)
01451     {
01452       register ssize_t
01453         i;
01454 
01455       if (GetPixelMask(image,p) != 0)
01456         {
01457           p+=GetPixelChannels(image);
01458           continue;
01459         }
01460       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
01461       {
01462         PixelChannel
01463           channel;
01464 
01465         PixelTrait
01466           traits;
01467 
01468         channel=GetPixelChannelMapChannel(image,i);
01469         traits=GetPixelChannelMapTraits(image,channel);
01470         if (traits == UndefinedPixelTrait)
01471           continue;
01472         if ((traits & UpdatePixelTrait) == 0)
01473           continue;
01474 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01475         #pragma omp critical (MagickCore_GetImageRange)
01476 #endif
01477         {
01478           if ((double) p[i] < *minima)
01479             *minima=(double) p[i];
01480           if ((double) p[i] > *maxima)
01481             *maxima=(double) p[i];
01482         }
01483       }
01484       p+=GetPixelChannels(image);
01485     }
01486   }
01487   image_view=DestroyCacheView(image_view);
01488   return(status);
01489 }
01490 
01491 /*
01492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01493 %                                                                             %
01494 %                                                                             %
01495 %                                                                             %
01496 %   G e t I m a g e S t a t i s t i c s                                       %
01497 %                                                                             %
01498 %                                                                             %
01499 %                                                                             %
01500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01501 %
01502 %  GetImageStatistics() returns statistics for each channel in the
01503 %  image.  The statistics include the channel depth, its minima, maxima, mean,
01504 %  standard deviation, kurtosis and skewness.  You can access the red channel
01505 %  mean, for example, like this:
01506 %
01507 %      channel_statistics=GetImageStatistics(image,exception);
01508 %      red_mean=channel_statistics[RedPixelChannel].mean;
01509 %
01510 %  Use MagickRelinquishMemory() to free the statistics buffer.
01511 %
01512 %  The format of the GetImageStatistics method is:
01513 %
01514 %      ChannelStatistics *GetImageStatistics(const Image *image,
01515 %        ExceptionInfo *exception)
01516 %
01517 %  A description of each parameter follows:
01518 %
01519 %    o image: the image.
01520 %
01521 %    o exception: return any errors or warnings in this structure.
01522 %
01523 */
01524 
01525 static size_t GetImageChannels(const Image *image)
01526 {
01527   register ssize_t
01528     i;
01529 
01530   size_t
01531     channels;
01532 
01533   channels=0;
01534   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
01535   {
01536     PixelChannel
01537       channel;
01538 
01539     PixelTrait
01540       traits;
01541 
01542     channel=GetPixelChannelMapChannel(image,i);
01543     traits=GetPixelChannelMapTraits(image,channel);
01544     if ((traits & UpdatePixelTrait) != 0)
01545       channels++;
01546   }
01547   return(channels);
01548 }
01549 
01550 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
01551   ExceptionInfo *exception)
01552 {
01553   ChannelStatistics
01554     *channel_statistics;
01555 
01556   double
01557     area;
01558 
01559   MagickStatusType
01560     status;
01561 
01562   QuantumAny
01563     range;
01564 
01565   register ssize_t
01566     i;
01567 
01568   size_t
01569     channels,
01570     depth;
01571 
01572   ssize_t
01573     y;
01574 
01575   assert(image != (Image *) NULL);
01576   assert(image->signature == MagickSignature);
01577   if (image->debug != MagickFalse)
01578     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01579   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
01580     MaxPixelChannels+1,sizeof(*channel_statistics));
01581   if (channel_statistics == (ChannelStatistics *) NULL)
01582     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
01583   (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
01584     sizeof(*channel_statistics));
01585   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
01586   {
01587     channel_statistics[i].depth=1;
01588     channel_statistics[i].maxima=(-MagickHuge);
01589     channel_statistics[i].minima=MagickHuge;
01590   }
01591   for (y=0; y < (ssize_t) image->rows; y++)
01592   {
01593     register const Quantum
01594       *restrict p;
01595 
01596     register ssize_t
01597       x;
01598 
01599     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
01600     if (p == (const Quantum *) NULL)
01601       break;
01602     for (x=0; x < (ssize_t) image->columns; x++)
01603     {
01604       register ssize_t
01605         i;
01606 
01607       if (GetPixelMask(image,p) != 0)
01608         {
01609           p+=GetPixelChannels(image);
01610           continue;
01611         }
01612       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
01613       {
01614         PixelChannel
01615           channel;
01616 
01617         PixelTrait
01618           traits;
01619 
01620         channel=GetPixelChannelMapChannel(image,i);
01621         traits=GetPixelChannelMapTraits(image,channel);
01622         if (traits == UndefinedPixelTrait)
01623           continue;
01624         if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
01625           {
01626             depth=channel_statistics[channel].depth;
01627             range=GetQuantumRange(depth);
01628             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
01629               range) ? MagickTrue : MagickFalse;
01630             if (status != MagickFalse)
01631               {
01632                 channel_statistics[channel].depth++;
01633                 continue;
01634               }
01635           }
01636         if ((double) p[i] < channel_statistics[channel].minima)
01637           channel_statistics[channel].minima=(double) p[i];
01638         if ((double) p[i] > channel_statistics[channel].maxima)
01639           channel_statistics[channel].maxima=(double) p[i];
01640         channel_statistics[channel].sum+=p[i];
01641         channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
01642         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
01643         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
01644           p[i];
01645       }
01646       p+=GetPixelChannels(image);
01647     }
01648   }
01649   area=(double) image->columns*image->rows;
01650   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
01651   {
01652     channel_statistics[i].sum/=area;
01653     channel_statistics[i].sum_squared/=area;
01654     channel_statistics[i].sum_cubed/=area;
01655     channel_statistics[i].sum_fourth_power/=area;
01656     channel_statistics[i].mean=channel_statistics[i].sum;
01657     channel_statistics[i].variance=channel_statistics[i].sum_squared;
01658     channel_statistics[i].standard_deviation=sqrt(
01659       channel_statistics[i].variance-(channel_statistics[i].mean*
01660       channel_statistics[i].mean));
01661   }
01662   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
01663   {
01664     channel_statistics[CompositePixelChannel].depth=(size_t) EvaluateMax(
01665       (double) channel_statistics[CompositePixelChannel].depth,(double)
01666       channel_statistics[i].depth);
01667     channel_statistics[CompositePixelChannel].minima=MagickMin(
01668       channel_statistics[CompositePixelChannel].minima,
01669       channel_statistics[i].minima);
01670     channel_statistics[CompositePixelChannel].maxima=EvaluateMax(
01671       channel_statistics[CompositePixelChannel].maxima,
01672       channel_statistics[i].maxima);
01673     channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
01674     channel_statistics[CompositePixelChannel].sum_squared+=
01675       channel_statistics[i].sum_squared;
01676     channel_statistics[CompositePixelChannel].sum_cubed+=
01677       channel_statistics[i].sum_cubed;
01678     channel_statistics[CompositePixelChannel].sum_fourth_power+=
01679       channel_statistics[i].sum_fourth_power;
01680     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
01681     channel_statistics[CompositePixelChannel].variance+=
01682       channel_statistics[i].variance-channel_statistics[i].mean*
01683       channel_statistics[i].mean;
01684     channel_statistics[CompositePixelChannel].standard_deviation+=
01685       channel_statistics[i].variance-channel_statistics[i].mean*
01686       channel_statistics[i].mean;
01687   }
01688   channels=GetImageChannels(image);
01689   channel_statistics[CompositePixelChannel].sum/=channels;
01690   channel_statistics[CompositePixelChannel].sum_squared/=channels;
01691   channel_statistics[CompositePixelChannel].sum_cubed/=channels;
01692   channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
01693   channel_statistics[CompositePixelChannel].mean/=channels;
01694   channel_statistics[CompositePixelChannel].variance/=channels;
01695   channel_statistics[CompositePixelChannel].standard_deviation=
01696     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
01697   channel_statistics[CompositePixelChannel].kurtosis/=channels;
01698   channel_statistics[CompositePixelChannel].skewness/=channels;
01699   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
01700   {
01701     if (channel_statistics[i].standard_deviation == 0.0)
01702       continue;
01703     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
01704       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
01705       channel_statistics[i].mean*channel_statistics[i].mean*
01706       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
01707       channel_statistics[i].standard_deviation*
01708       channel_statistics[i].standard_deviation);
01709     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
01710       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
01711       channel_statistics[i].mean*channel_statistics[i].mean*
01712       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
01713       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
01714       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
01715       channel_statistics[i].standard_deviation*
01716       channel_statistics[i].standard_deviation*
01717       channel_statistics[i].standard_deviation)-3.0;
01718   }
01719   return(channel_statistics);
01720 }
01721 
01722 /*
01723 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01724 %                                                                             %
01725 %                                                                             %
01726 %                                                                             %
01727 %     S t a t i s t i c I m a g e                                             %
01728 %                                                                             %
01729 %                                                                             %
01730 %                                                                             %
01731 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01732 %
01733 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
01734 %  the neighborhood of the specified width and height.
01735 %
01736 %  The format of the StatisticImage method is:
01737 %
01738 %      Image *StatisticImage(const Image *image,const StatisticType type,
01739 %        const size_t width,const size_t height,ExceptionInfo *exception)
01740 %
01741 %  A description of each parameter follows:
01742 %
01743 %    o image: the image.
01744 %
01745 %    o type: the statistic type (median, mode, etc.).
01746 %
01747 %    o width: the width of the pixel neighborhood.
01748 %
01749 %    o height: the height of the pixel neighborhood.
01750 %
01751 %    o exception: return any errors or warnings in this structure.
01752 %
01753 */
01754 
01755 typedef struct _SkipNode
01756 {
01757   size_t
01758     next[9],
01759     count,
01760     signature;
01761 } SkipNode;
01762 
01763 typedef struct _SkipList
01764 {
01765   ssize_t
01766     level;
01767 
01768   SkipNode
01769     *nodes;
01770 } SkipList;
01771 
01772 typedef struct _PixelList
01773 {
01774   size_t
01775     length,
01776     seed;
01777 
01778   SkipList
01779     skip_list;
01780 
01781   size_t
01782     signature;
01783 } PixelList;
01784 
01785 static PixelList *DestroyPixelList(PixelList *pixel_list)
01786 {
01787   if (pixel_list == (PixelList *) NULL)
01788     return((PixelList *) NULL);
01789   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
01790     pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
01791       pixel_list->skip_list.nodes);
01792   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
01793   return(pixel_list);
01794 }
01795 
01796 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
01797 {
01798   register ssize_t
01799     i;
01800 
01801   assert(pixel_list != (PixelList **) NULL);
01802   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
01803     if (pixel_list[i] != (PixelList *) NULL)
01804       pixel_list[i]=DestroyPixelList(pixel_list[i]);
01805   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
01806   return(pixel_list);
01807 }
01808 
01809 static PixelList *AcquirePixelList(const size_t width,const size_t height)
01810 {
01811   PixelList
01812     *pixel_list;
01813 
01814   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
01815   if (pixel_list == (PixelList *) NULL)
01816     return(pixel_list);
01817   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
01818   pixel_list->length=width*height;
01819   pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
01820     sizeof(*pixel_list->skip_list.nodes));
01821   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
01822     return(DestroyPixelList(pixel_list));
01823   (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
01824     sizeof(*pixel_list->skip_list.nodes));
01825   pixel_list->signature=MagickSignature;
01826   return(pixel_list);
01827 }
01828 
01829 static PixelList **AcquirePixelListThreadSet(const size_t width,
01830   const size_t height)
01831 {
01832   PixelList
01833     **pixel_list;
01834 
01835   register ssize_t
01836     i;
01837 
01838   size_t
01839     number_threads;
01840 
01841   number_threads=GetOpenMPMaximumThreads();
01842   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
01843     sizeof(*pixel_list));
01844   if (pixel_list == (PixelList **) NULL)
01845     return((PixelList **) NULL);
01846   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
01847   for (i=0; i < (ssize_t) number_threads; i++)
01848   {
01849     pixel_list[i]=AcquirePixelList(width,height);
01850     if (pixel_list[i] == (PixelList *) NULL)
01851       return(DestroyPixelListThreadSet(pixel_list));
01852   }
01853   return(pixel_list);
01854 }
01855 
01856 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
01857 {
01858   register SkipList
01859     *p;
01860 
01861   register ssize_t
01862     level;
01863 
01864   size_t
01865     search,
01866     update[9];
01867 
01868   /*
01869     Initialize the node.
01870   */
01871   p=(&pixel_list->skip_list);
01872   p->nodes[color].signature=pixel_list->signature;
01873   p->nodes[color].count=1;
01874   /*
01875     Determine where it belongs in the list.
01876   */
01877   search=65536UL;
01878   for (level=p->level; level >= 0; level--)
01879   {
01880     while (p->nodes[search].next[level] < color)
01881       search=p->nodes[search].next[level];
01882     update[level]=search;
01883   }
01884   /*
01885     Generate a pseudo-random level for this node.
01886   */
01887   for (level=0; ; level++)
01888   {
01889     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
01890     if ((pixel_list->seed & 0x300) != 0x300)
01891       break;
01892   }
01893   if (level > 8)
01894     level=8;
01895   if (level > (p->level+2))
01896     level=p->level+2;
01897   /*
01898     If we're raising the list's level, link back to the root node.
01899   */
01900   while (level > p->level)
01901   {
01902     p->level++;
01903     update[p->level]=65536UL;
01904   }
01905   /*
01906     Link the node into the skip-list.
01907   */
01908   do
01909   {
01910     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
01911     p->nodes[update[level]].next[level]=color;
01912   } while (level-- > 0);
01913 }
01914 
01915 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
01916 {
01917   register SkipList
01918     *p;
01919 
01920   size_t
01921     color,
01922     maximum;
01923 
01924   ssize_t
01925     count;
01926 
01927   /*
01928     Find the maximum value for each of the color.
01929   */
01930   p=(&pixel_list->skip_list);
01931   color=65536L;
01932   count=0;
01933   maximum=p->nodes[color].next[0];
01934   do
01935   {
01936     color=p->nodes[color].next[0];
01937     if (color > maximum)
01938       maximum=color;
01939     count+=p->nodes[color].count;
01940   } while (count < (ssize_t) pixel_list->length);
01941   *pixel=ScaleShortToQuantum((unsigned short) maximum);
01942 }
01943 
01944 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
01945 {
01946   MagickRealType
01947     sum;
01948 
01949   register SkipList
01950     *p;
01951 
01952   size_t
01953     color;
01954 
01955   ssize_t
01956     count;
01957 
01958   /*
01959     Find the mean value for each of the color.
01960   */
01961   p=(&pixel_list->skip_list);
01962   color=65536L;
01963   count=0;
01964   sum=0.0;
01965   do
01966   {
01967     color=p->nodes[color].next[0];
01968     sum+=(MagickRealType) p->nodes[color].count*color;
01969     count+=p->nodes[color].count;
01970   } while (count < (ssize_t) pixel_list->length);
01971   sum/=pixel_list->length;
01972   *pixel=ScaleShortToQuantum((unsigned short) sum);
01973 }
01974 
01975 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
01976 {
01977   register SkipList
01978     *p;
01979 
01980   size_t
01981     color;
01982 
01983   ssize_t
01984     count;
01985 
01986   /*
01987     Find the median value for each of the color.
01988   */
01989   p=(&pixel_list->skip_list);
01990   color=65536L;
01991   count=0;
01992   do
01993   {
01994     color=p->nodes[color].next[0];
01995     count+=p->nodes[color].count;
01996   } while (count <= (ssize_t) (pixel_list->length >> 1));
01997   *pixel=ScaleShortToQuantum((unsigned short) color);
01998 }
01999 
02000 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
02001 {
02002   register SkipList
02003     *p;
02004 
02005   size_t
02006     color,
02007     minimum;
02008 
02009   ssize_t
02010     count;
02011 
02012   /*
02013     Find the minimum value for each of the color.
02014   */
02015   p=(&pixel_list->skip_list);
02016   count=0;
02017   color=65536UL;
02018   minimum=p->nodes[color].next[0];
02019   do
02020   {
02021     color=p->nodes[color].next[0];
02022     if (color < minimum)
02023       minimum=color;
02024     count+=p->nodes[color].count;
02025   } while (count < (ssize_t) pixel_list->length);
02026   *pixel=ScaleShortToQuantum((unsigned short) minimum);
02027 }
02028 
02029 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
02030 {
02031   register SkipList
02032     *p;
02033 
02034   size_t
02035     color,
02036     max_count,
02037     mode;
02038 
02039   ssize_t
02040     count;
02041 
02042   /*
02043     Make each pixel the 'predominant color' of the specified neighborhood.
02044   */
02045   p=(&pixel_list->skip_list);
02046   color=65536L;
02047   mode=color;
02048   max_count=p->nodes[mode].count;
02049   count=0;
02050   do
02051   {
02052     color=p->nodes[color].next[0];
02053     if (p->nodes[color].count > max_count)
02054       {
02055         mode=color;
02056         max_count=p->nodes[mode].count;
02057       }
02058     count+=p->nodes[color].count;
02059   } while (count < (ssize_t) pixel_list->length);
02060   *pixel=ScaleShortToQuantum((unsigned short) mode);
02061 }
02062 
02063 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
02064 {
02065   register SkipList
02066     *p;
02067 
02068   size_t
02069     color,
02070     next,
02071     previous;
02072 
02073   ssize_t
02074     count;
02075 
02076   /*
02077     Finds the non peak value for each of the colors.
02078   */
02079   p=(&pixel_list->skip_list);
02080   color=65536L;
02081   next=p->nodes[color].next[0];
02082   count=0;
02083   do
02084   {
02085     previous=color;
02086     color=next;
02087     next=p->nodes[color].next[0];
02088     count+=p->nodes[color].count;
02089   } while (count <= (ssize_t) (pixel_list->length >> 1));
02090   if ((previous == 65536UL) && (next != 65536UL))
02091     color=next;
02092   else
02093     if ((previous != 65536UL) && (next == 65536UL))
02094       color=previous;
02095   *pixel=ScaleShortToQuantum((unsigned short) color);
02096 }
02097 
02098 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
02099   Quantum *pixel)
02100 {
02101   MagickRealType
02102     sum,
02103     sum_squared;
02104 
02105   register SkipList
02106     *p;
02107 
02108   size_t
02109     color;
02110 
02111   ssize_t
02112     count;
02113 
02114   /*
02115     Find the standard-deviation value for each of the color.
02116   */
02117   p=(&pixel_list->skip_list);
02118   color=65536L;
02119   count=0;
02120   sum=0.0;
02121   sum_squared=0.0;
02122   do
02123   {
02124     register ssize_t
02125       i;
02126 
02127     color=p->nodes[color].next[0];
02128     sum+=(MagickRealType) p->nodes[color].count*color;
02129     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
02130       sum_squared+=((MagickRealType) color)*((MagickRealType) color);
02131     count+=p->nodes[color].count;
02132   } while (count < (ssize_t) pixel_list->length);
02133   sum/=pixel_list->length;
02134   sum_squared/=pixel_list->length;
02135   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
02136 }
02137 
02138 static inline void InsertPixelList(const Image *image,const Quantum pixel,
02139   PixelList *pixel_list)
02140 {
02141   size_t
02142     signature;
02143 
02144   unsigned short
02145     index;
02146 
02147   index=ScaleQuantumToShort(pixel);
02148   signature=pixel_list->skip_list.nodes[index].signature;
02149   if (signature == pixel_list->signature)
02150     {
02151       pixel_list->skip_list.nodes[index].count++;
02152       return;
02153     }
02154   AddNodePixelList(pixel_list,index);
02155 }
02156 
02157 static inline MagickRealType MagickAbsoluteValue(const MagickRealType x)
02158 {
02159   if (x < 0)
02160     return(-x);
02161   return(x);
02162 }
02163 
02164 static inline size_t MagickMax(const size_t x,const size_t y)
02165 {
02166   if (x > y)
02167     return(x);
02168   return(y);
02169 }
02170 
02171 static void ResetPixelList(PixelList *pixel_list)
02172 {
02173   int
02174     level;
02175 
02176   register SkipNode
02177     *root;
02178 
02179   register SkipList
02180     *p;
02181 
02182   /*
02183     Reset the skip-list.
02184   */
02185   p=(&pixel_list->skip_list);
02186   root=p->nodes+65536UL;
02187   p->level=0;
02188   for (level=0; level < 9; level++)
02189     root->next[level]=65536UL;
02190   pixel_list->seed=pixel_list->signature++;
02191 }
02192 
02193 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
02194   const size_t width,const size_t height,ExceptionInfo *exception)
02195 {
02196 #define StatisticImageTag  "Statistic/Image"
02197 
02198   CacheView
02199     *image_view,
02200     *statistic_view;
02201 
02202   Image
02203     *statistic_image;
02204 
02205   MagickBooleanType
02206     status;
02207 
02208   MagickOffsetType
02209     progress;
02210 
02211   PixelList
02212     **restrict pixel_list;
02213 
02214   ssize_t
02215     center,
02216     y;
02217 
02218   /*
02219     Initialize statistics image attributes.
02220   */
02221   assert(image != (Image *) NULL);
02222   assert(image->signature == MagickSignature);
02223   if (image->debug != MagickFalse)
02224     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
02225   assert(exception != (ExceptionInfo *) NULL);
02226   assert(exception->signature == MagickSignature);
02227   statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
02228     exception);
02229   if (statistic_image == (Image *) NULL)
02230     return((Image *) NULL);
02231   status=SetImageStorageClass(statistic_image,DirectClass,exception);
02232   if (status == MagickFalse)
02233     {
02234       statistic_image=DestroyImage(statistic_image);
02235       return((Image *) NULL);
02236     }
02237   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
02238   if (pixel_list == (PixelList **) NULL)
02239     {
02240       statistic_image=DestroyImage(statistic_image);
02241       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
02242     }
02243   /*
02244     Make each pixel the min / max / median / mode / etc. of the neighborhood.
02245   */
02246   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
02247     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
02248   status=MagickTrue;
02249   progress=0;
02250   image_view=AcquireCacheView(image);
02251   statistic_view=AcquireCacheView(statistic_image);
02252 #if defined(MAGICKCORE_OPENMP_SUPPORT)
02253   #pragma omp parallel for schedule(static,4) shared(progress,status)
02254 #endif
02255   for (y=0; y < (ssize_t) statistic_image->rows; y++)
02256   {
02257     const int
02258       id = GetOpenMPThreadId();
02259 
02260     register const Quantum
02261       *restrict p;
02262 
02263     register Quantum
02264       *restrict q;
02265 
02266     register ssize_t
02267       x;
02268 
02269     if (status == MagickFalse)
02270       continue;
02271     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
02272       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
02273       MagickMax(height,1),exception);
02274     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
02275     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
02276       {
02277         status=MagickFalse;
02278         continue;
02279       }
02280     for (x=0; x < (ssize_t) statistic_image->columns; x++)
02281     {
02282       register ssize_t
02283         i;
02284 
02285       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
02286       {
02287         PixelChannel
02288           channel;
02289 
02290         PixelTrait
02291           statistic_traits,
02292           traits;
02293 
02294         Quantum
02295           pixel;
02296 
02297         register const Quantum
02298           *restrict pixels;
02299 
02300         register ssize_t
02301           u;
02302 
02303         ssize_t
02304           v;
02305 
02306         channel=GetPixelChannelMapChannel(image,i);
02307         traits=GetPixelChannelMapTraits(image,channel);
02308         statistic_traits=GetPixelChannelMapTraits(statistic_image,channel);
02309         if ((traits == UndefinedPixelTrait) ||
02310             (statistic_traits == UndefinedPixelTrait))
02311           continue;
02312         if (((statistic_traits & CopyPixelTrait) != 0) ||
02313             (GetPixelMask(image,p) != 0))
02314           {
02315             SetPixelChannel(statistic_image,channel,p[center+i],q);
02316             continue;
02317           }
02318         pixels=p;
02319         ResetPixelList(pixel_list[id]);
02320         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
02321         {
02322           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
02323           {
02324             InsertPixelList(image,pixels[i],pixel_list[id]);
02325             pixels+=GetPixelChannels(image);
02326           }
02327           pixels+=image->columns*GetPixelChannels(image);
02328         }
02329         switch (type)
02330         {
02331           case GradientStatistic:
02332           {
02333             MagickRealType
02334               maximum,
02335               minimum;
02336 
02337             GetMinimumPixelList(pixel_list[id],&pixel);
02338             minimum=(MagickRealType) pixel;
02339             GetMaximumPixelList(pixel_list[id],&pixel);
02340             maximum=(MagickRealType) pixel;
02341             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
02342             break;
02343           }
02344           case MaximumStatistic:
02345           {
02346             GetMaximumPixelList(pixel_list[id],&pixel);
02347             break;
02348           }
02349           case MeanStatistic:
02350           {
02351             GetMeanPixelList(pixel_list[id],&pixel);
02352             break;
02353           }
02354           case MedianStatistic:
02355           default:
02356           {
02357             GetMedianPixelList(pixel_list[id],&pixel);
02358             break;
02359           }
02360           case MinimumStatistic:
02361           {
02362             GetMinimumPixelList(pixel_list[id],&pixel);
02363             break;
02364           }
02365           case ModeStatistic:
02366           {
02367             GetModePixelList(pixel_list[id],&pixel);
02368             break;
02369           }
02370           case NonpeakStatistic:
02371           {
02372             GetNonpeakPixelList(pixel_list[id],&pixel);
02373             break;
02374           }
02375           case StandardDeviationStatistic:
02376           {
02377             GetStandardDeviationPixelList(pixel_list[id],&pixel);
02378             break;
02379           }
02380         }
02381         SetPixelChannel(statistic_image,channel,pixel,q);
02382       }
02383       p+=GetPixelChannels(image);
02384       q+=GetPixelChannels(statistic_image);
02385     }
02386     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
02387       status=MagickFalse;
02388     if (image->progress_monitor != (MagickProgressMonitor) NULL)
02389       {
02390         MagickBooleanType
02391           proceed;
02392 
02393 #if defined(MAGICKCORE_OPENMP_SUPPORT)
02394   #pragma omp critical (MagickCore_StatisticImage)
02395 #endif
02396         proceed=SetImageProgress(image,StatisticImageTag,progress++,
02397           image->rows);
02398         if (proceed == MagickFalse)
02399           status=MagickFalse;
02400       }
02401   }
02402   statistic_view=DestroyCacheView(statistic_view);
02403   image_view=DestroyCacheView(image_view);
02404   pixel_list=DestroyPixelListThreadSet(pixel_list);
02405   return(statistic_image);
02406 }