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 Methods                           %
00014 %                                                                             %
00015 %                              Software Design                                %
00016 %                                John Cristy                                  %
00017 %                                 July 1992                                   %
00018 %                                                                             %
00019 %                                                                             %
00020 %  Copyright 1999-2009 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 "magick/studio.h"
00044 #include "magick/property.h"
00045 #include "magick/animate.h"
00046 #include "magick/blob.h"
00047 #include "magick/blob-private.h"
00048 #include "magick/cache.h"
00049 #include "magick/cache-private.h"
00050 #include "magick/cache-view.h"
00051 #include "magick/client.h"
00052 #include "magick/color.h"
00053 #include "magick/color-private.h"
00054 #include "magick/colorspace.h"
00055 #include "magick/colorspace-private.h"
00056 #include "magick/composite.h"
00057 #include "magick/composite-private.h"
00058 #include "magick/compress.h"
00059 #include "magick/constitute.h"
00060 #include "magick/deprecate.h"
00061 #include "magick/display.h"
00062 #include "magick/draw.h"
00063 #include "magick/enhance.h"
00064 #include "magick/exception.h"
00065 #include "magick/exception-private.h"
00066 #include "magick/gem.h"
00067 #include "magick/geometry.h"
00068 #include "magick/list.h"
00069 #include "magick/image-private.h"
00070 #include "magick/magic.h"
00071 #include "magick/magick.h"
00072 #include "magick/memory_.h"
00073 #include "magick/module.h"
00074 #include "magick/monitor.h"
00075 #include "magick/monitor-private.h"
00076 #include "magick/option.h"
00077 #include "magick/paint.h"
00078 #include "magick/pixel-private.h"
00079 #include "magick/profile.h"
00080 #include "magick/quantize.h"
00081 #include "magick/random_.h"
00082 #include "magick/segment.h"
00083 #include "magick/semaphore.h"
00084 #include "magick/signature-private.h"
00085 #include "magick/statistic.h"
00086 #include "magick/string_.h"
00087 #include "magick/thread-private.h"
00088 #include "magick/timer.h"
00089 #include "magick/utility.h"
00090 #include "magick/version.h"
00091 
00092 /*
00093 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00094 %                                                                             %
00095 %                                                                             %
00096 %                                                                             %
00097 %     A v e r a g e I m a g e s                                               %
00098 %                                                                             %
00099 %                                                                             %
00100 %                                                                             %
00101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00102 %
00103 %  AverageImages() takes a set of images and averages them together.  Each
00104 %  image in the set must have the same width and height.  AverageImages()
00105 %  returns a single image with each corresponding pixel component of each
00106 %  image averaged.   On failure, a NULL image is returned and exception
00107 %  describes the reason for the failure.
00108 %
00109 %  The format of the AverageImages method is:
00110 %
00111 %      Image *AverageImages(Image *image,ExceptionInfo *exception)
00112 %
00113 %  A description of each parameter follows:
00114 %
00115 %    o image: the image sequence.
00116 %
00117 %    o exception: return any errors or warnings in this structure.
00118 %
00119 */
00120 
00121 static MagickPixelPacket **DestroyPixelThreadSet(MagickPixelPacket **pixels)
00122 {
00123   register long
00124     i;
00125 
00126   assert(pixels != (MagickPixelPacket **) NULL);
00127   for (i=0; i < (long) GetOpenMPMaximumThreads(); i++)
00128     if (pixels[i] != (MagickPixelPacket *) NULL)
00129       pixels[i]=(MagickPixelPacket *) RelinquishMagickMemory(pixels[i]);
00130   pixels=(MagickPixelPacket **) RelinquishAlignedMemory(pixels);
00131   return(pixels);
00132 }
00133 
00134 static MagickPixelPacket **AcquirePixelThreadSet(const Image *image)
00135 {
00136   register long
00137     i,
00138     j;
00139 
00140   MagickPixelPacket
00141     **pixels;
00142 
00143   unsigned long
00144     number_threads;
00145 
00146   number_threads=GetOpenMPMaximumThreads();
00147   pixels=(MagickPixelPacket **) AcquireAlignedMemory(number_threads,
00148     sizeof(*pixels));
00149   if (pixels == (MagickPixelPacket **) NULL)
00150     return((MagickPixelPacket **) NULL);
00151   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
00152   for (i=0; i < (long) number_threads; i++)
00153   {
00154     pixels[i]=(MagickPixelPacket *) AcquireQuantumMemory(image->columns,
00155       sizeof(**pixels));
00156     if (pixels[i] == (MagickPixelPacket *) NULL)
00157       return(DestroyPixelThreadSet(pixels));
00158     for (j=0; j < (long) image->columns; j++)
00159       GetMagickPixelPacket(image,&pixels[i][j]);
00160   }
00161   return(pixels);
00162 }
00163 
00164 MagickExport Image *AverageImages(const Image *image,ExceptionInfo *exception)
00165 {
00166 #define AverageImageTag  "Average/Image"
00167 
00168   CacheView
00169     *average_view;
00170 
00171   const Image
00172     *next;
00173 
00174   Image
00175     *average_image;
00176 
00177   long
00178     progress,
00179     y;
00180 
00181   MagickBooleanType
00182     status;
00183 
00184   MagickPixelPacket
00185     **average_pixels,
00186     zero;
00187 
00188   unsigned long
00189     number_images;
00190 
00191   /*
00192     Ensure the image are the same size.
00193   */
00194   assert(image != (Image *) NULL);
00195   assert(image->signature == MagickSignature);
00196   if (image->debug != MagickFalse)
00197     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00198   assert(exception != (ExceptionInfo *) NULL);
00199   assert(exception->signature == MagickSignature);
00200   for (next=image; next != (Image *) NULL; next=GetNextImageInList(next))
00201     if ((next->columns != image->columns) || (next->rows != image->rows))
00202       ThrowImageException(OptionError,"ImageWidthsOrHeightsDiffer");
00203   /*
00204     Initialize average next attributes.
00205   */
00206   average_image=CloneImage(image,image->columns,image->rows,MagickTrue,
00207     exception);
00208   if (average_image == (Image *) NULL)
00209     return((Image *) NULL);
00210   if (SetImageStorageClass(average_image,DirectClass) == MagickFalse)
00211     {
00212       InheritException(exception,&average_image->exception);
00213       average_image=DestroyImage(average_image);
00214       return((Image *) NULL);
00215     }
00216   average_pixels=AcquirePixelThreadSet(image);
00217   if (average_pixels == (MagickPixelPacket **) NULL)
00218     {
00219       average_image=DestroyImage(average_image);
00220       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
00221     }
00222   /*
00223     Average image pixels.
00224   */
00225   status=MagickTrue;
00226   progress=0;
00227   GetMagickPixelPacket(image,&zero);
00228   number_images=GetImageListLength(image);
00229   average_view=AcquireCacheView(average_image);
00230 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00231   #pragma omp parallel for schedule(dynamic) shared(progress,status)
00232 #endif
00233   for (y=0; y < (long) average_image->rows; y++)
00234   {
00235     CacheView
00236       *image_view;
00237 
00238     const Image
00239       *next;
00240 
00241     MagickPixelPacket
00242       pixel;
00243 
00244     register IndexPacket
00245       *__restrict average_indexes;
00246 
00247     register long
00248       i,
00249       id,
00250       x;
00251 
00252     register MagickPixelPacket
00253       *average_pixel;
00254 
00255     register PixelPacket
00256       *__restrict q;
00257 
00258     if (status == MagickFalse)
00259       continue;
00260     q=QueueCacheViewAuthenticPixels(average_view,0,y,average_image->columns,1,
00261       exception);
00262     if (q == (PixelPacket *) NULL)
00263       {
00264         status=MagickFalse;
00265         continue;
00266       }
00267     average_indexes=GetCacheViewAuthenticIndexQueue(average_view);
00268     pixel=zero;
00269     id=GetOpenMPThreadId();
00270     average_pixel=average_pixels[id];
00271     for (x=0; x < (long) average_image->columns; x++)
00272       average_pixel[x]=zero;
00273     next=image;
00274     for (i=0; i < (long) number_images; i++)
00275     {
00276       register const IndexPacket
00277         *indexes;
00278 
00279       register const PixelPacket
00280         *p;
00281 
00282       image_view=AcquireCacheView(next);
00283       p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
00284       if (p == (const PixelPacket *) NULL)
00285         {
00286           image_view=DestroyCacheView(image_view);
00287           break;
00288         }
00289       indexes=GetCacheViewVirtualIndexQueue(image_view);
00290       for (x=0; x < (long) next->columns; x++)
00291       {
00292         SetMagickPixelPacket(next,p,indexes+x,&pixel);
00293         average_pixel[x].red+=QuantumScale*pixel.red;
00294         average_pixel[x].green+=QuantumScale*pixel.green;
00295         average_pixel[x].blue+=QuantumScale*pixel.blue;
00296         average_pixel[x].opacity+=QuantumScale*pixel.opacity;
00297         if (average_image->colorspace == CMYKColorspace)
00298           average_pixel[x].index+=QuantumScale*pixel.index;
00299         p++;
00300       }
00301       image_view=DestroyCacheView(image_view);
00302       next=GetNextImageInList(next);
00303     }
00304     for (x=0; x < (long) average_image->columns; x++)
00305     {
00306       average_pixel[x].red=(MagickRealType) (QuantumRange*
00307         average_pixel[x].red/number_images);
00308       average_pixel[x].green=(MagickRealType) (QuantumRange*
00309         average_pixel[x].green/number_images);
00310       average_pixel[x].blue=(MagickRealType) (QuantumRange*
00311         average_pixel[x].blue/number_images);
00312       average_pixel[x].opacity=(MagickRealType) (QuantumRange*
00313         average_pixel[x].opacity/number_images);
00314       if (average_image->colorspace == CMYKColorspace)
00315         average_pixel[x].index=(MagickRealType) (QuantumRange*
00316           average_pixel[x].index/number_images);
00317       SetPixelPacket(average_image,&average_pixel[x],q,average_indexes+x);
00318       q++;
00319     }
00320     if (SyncCacheViewAuthenticPixels(average_view,exception) == MagickFalse)
00321       status=MagickFalse;
00322     if (image->progress_monitor != (MagickProgressMonitor) NULL)
00323       {
00324         MagickBooleanType
00325           proceed;
00326 
00327 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00328         #pragma omp critical (MagickCore_AverageImages)
00329 #endif
00330         proceed=SetImageProgress(image,AverageImageTag,progress++,
00331           average_image->rows);
00332         if (proceed == MagickFalse)
00333           status=MagickFalse;
00334       }
00335   }
00336   average_view=DestroyCacheView(average_view);
00337   average_pixels=DestroyPixelThreadSet(average_pixels);
00338   if (status == MagickFalse)
00339     average_image=DestroyImage(average_image);
00340   return(average_image);
00341 }
00342 
00343 /*
00344 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00345 %                                                                             %
00346 %                                                                             %
00347 %                                                                             %
00348 +   G e t I m a g e C h a n n e l E x t r e m a                               %
00349 %                                                                             %
00350 %                                                                             %
00351 %                                                                             %
00352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00353 %
00354 %  GetImageChannelExtrema() returns the extrema of one or more image channels.
00355 %
00356 %  The format of the GetImageChannelExtrema method is:
00357 %
00358 %      MagickBooleanType GetImageChannelExtrema(const Image *image,
00359 %        const ChannelType channel,unsigned long *minima,unsigned long *maxima,
00360 %        ExceptionInfo *exception)
00361 %
00362 %  A description of each parameter follows:
00363 %
00364 %    o image: the image.
00365 %
00366 %    o channel: the channel.
00367 %
00368 %    o minima: the minimum value in the channel.
00369 %
00370 %    o maxima: the maximum value in the channel.
00371 %
00372 %    o exception: return any errors or warnings in this structure.
00373 %
00374 */
00375 
00376 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
00377   unsigned long *minima,unsigned long *maxima,ExceptionInfo *exception)
00378 {
00379   return(GetImageChannelExtrema(image,AllChannels,minima,maxima,exception));
00380 }
00381 
00382 MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
00383   const ChannelType channel,unsigned long *minima,unsigned long *maxima,
00384   ExceptionInfo *exception)
00385 {
00386   double
00387     max,
00388     min;
00389 
00390   MagickBooleanType
00391     status;
00392 
00393   assert(image != (Image *) NULL);
00394   assert(image->signature == MagickSignature);
00395   if (image->debug != MagickFalse)
00396     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00397   status=GetImageChannelRange(image,channel,&min,&max,exception);
00398   *minima=(unsigned long) (min+0.5);
00399   *maxima=(unsigned long) (max+0.5);
00400   return(status);
00401 }
00402 
00403 /*
00404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00405 %                                                                             %
00406 %                                                                             %
00407 %                                                                             %
00408 %   G e t I m a g e C h a n n e l M e a n                                     %
00409 %                                                                             %
00410 %                                                                             %
00411 %                                                                             %
00412 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00413 %
00414 %  GetImageChannelMean() returns the mean and standard deviation of one or more
00415 %  image channels.
00416 %
00417 %  The format of the GetImageChannelMean method is:
00418 %
00419 %      MagickBooleanType GetImageChannelMean(const Image *image,
00420 %        const ChannelType channel,double *mean,double *standard_deviation,
00421 %        ExceptionInfo *exception)
00422 %
00423 %  A description of each parameter follows:
00424 %
00425 %    o image: the image.
00426 %
00427 %    o channel: the channel.
00428 %
00429 %    o mean: the average value in the channel.
00430 %
00431 %    o standard_deviation: the standard deviation of the channel.
00432 %
00433 %    o exception: return any errors or warnings in this structure.
00434 %
00435 */
00436 
00437 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
00438   double *standard_deviation,ExceptionInfo *exception)
00439 {
00440   MagickBooleanType
00441     status;
00442 
00443   status=GetImageChannelMean(image,AllChannels,mean,standard_deviation,
00444     exception);
00445   return(status);
00446 }
00447 
00448 MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
00449   const ChannelType channel,double *mean,double *standard_deviation,
00450   ExceptionInfo *exception)
00451 {
00452   double
00453     area;
00454 
00455   long
00456     y;
00457 
00458   assert(image != (Image *) NULL);
00459   assert(image->signature == MagickSignature);
00460   if (image->debug != MagickFalse)
00461     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00462   *mean=0.0;
00463   *standard_deviation=0.0;
00464   area=0.0;
00465   for (y=0; y < (long) image->rows; y++)
00466   {
00467     register const IndexPacket
00468       *__restrict indexes;
00469 
00470     register const PixelPacket
00471       *__restrict p;
00472 
00473     register long
00474       x;
00475 
00476     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
00477     if (p == (const PixelPacket *) NULL)
00478       break;
00479     indexes=GetVirtualIndexQueue(image);
00480     for (x=0; x < (long) image->columns; x++)
00481     {
00482       if ((channel & RedChannel) != 0)
00483         {
00484           *mean+=p->red;
00485           *standard_deviation+=(double) p->red*p->red;
00486           area++;
00487         }
00488       if ((channel & GreenChannel) != 0)
00489         {
00490           *mean+=p->green;
00491           *standard_deviation+=(double) p->green*p->green;
00492           area++;
00493         }
00494       if ((channel & BlueChannel) != 0)
00495         {
00496           *mean+=p->blue;
00497           *standard_deviation+=(double) p->blue*p->blue;
00498           area++;
00499         }
00500       if ((channel & OpacityChannel) != 0)
00501         {
00502           *mean+=p->opacity;
00503           *standard_deviation+=(double) p->opacity*p->opacity;
00504           area++;
00505         }
00506       if (((channel & IndexChannel) != 0) &&
00507           (image->colorspace == CMYKColorspace))
00508         {
00509           *mean+=indexes[x];
00510           *standard_deviation+=(double) indexes[x]*indexes[x];
00511           area++;
00512         }
00513       p++;
00514     }
00515   }
00516   if (y < (long) image->rows)
00517     return(MagickFalse);
00518   if (area != 0)
00519     {
00520       *mean/=area;
00521       *standard_deviation/=area;
00522     }
00523   *standard_deviation=sqrt(*standard_deviation-(*mean*(*mean)));
00524   return(y == (long) image->rows ? MagickTrue : MagickFalse);
00525 }
00526 
00527 /*
00528 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00529 %                                                                             %
00530 %                                                                             %
00531 %                                                                             %
00532 %   G e t I m a g e C h a n n e l K u r t o s i s                             %
00533 %                                                                             %
00534 %                                                                             %
00535 %                                                                             %
00536 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00537 %
00538 %  GetImageChannelKurtosis() returns the kurtosis and skewness of one or more
00539 %  image channels.
00540 %
00541 %  The format of the GetImageChannelKurtosis method is:
00542 %
00543 %      MagickBooleanType GetImageChannelKurtosis(const Image *image,
00544 %        const ChannelType channel,double *kurtosis,double *skewness,
00545 %        ExceptionInfo *exception)
00546 %
00547 %  A description of each parameter follows:
00548 %
00549 %    o image: the image.
00550 %
00551 %    o channel: the channel.
00552 %
00553 %    o kurtosis: the kurtosis of the channel.
00554 %
00555 %    o skewness: the skewness of the channel.
00556 %
00557 %    o exception: return any errors or warnings in this structure.
00558 %
00559 */
00560 
00561 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
00562   double *kurtosis,double *skewness,ExceptionInfo *exception)
00563 {
00564   MagickBooleanType
00565     status;
00566 
00567   status=GetImageChannelKurtosis(image,AllChannels,kurtosis,skewness,
00568     exception);
00569   return(status);
00570 }
00571 
00572 MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
00573   const ChannelType channel,double *kurtosis,double *skewness,
00574   ExceptionInfo *exception)
00575 {
00576   double
00577     area,
00578     mean,
00579     standard_deviation,
00580     sum_squares,
00581     sum_cubes,
00582     sum_fourth_power;
00583 
00584   long
00585     y;
00586 
00587   assert(image != (Image *) NULL);
00588   assert(image->signature == MagickSignature);
00589   if (image->debug != MagickFalse)
00590     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00591   *kurtosis=0.0;
00592   *skewness=0.0;
00593   area=0.0;
00594   mean=0.0;
00595   standard_deviation=0.0;
00596   sum_squares=0.0;
00597   sum_cubes=0.0;
00598   sum_fourth_power=0.0;
00599   for (y=0; y < (long) image->rows; y++)
00600   {
00601     register const IndexPacket
00602       *__restrict indexes;
00603 
00604     register const PixelPacket
00605       *__restrict p;
00606 
00607     register long
00608       x;
00609 
00610     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
00611     if (p == (const PixelPacket *) NULL)
00612       break;
00613     indexes=GetVirtualIndexQueue(image);
00614     for (x=0; x < (long) image->columns; x++)
00615     {
00616       if ((channel & RedChannel) != 0)
00617         {
00618           mean+=p->red;
00619           sum_squares+=(double) p->red*p->red;
00620           sum_cubes+=(double) p->red*p->red*p->red;
00621           sum_fourth_power+=(double) p->red*p->red*p->red*p->red;
00622           area++;
00623         }
00624       if ((channel & GreenChannel) != 0)
00625         {
00626           mean+=p->green;
00627           sum_squares+=(double) p->green*p->green;
00628           sum_cubes+=(double) p->green*p->green*p->green;
00629           sum_fourth_power+=(double) p->green*p->green*p->green*p->green;
00630           area++;
00631         }
00632       if ((channel & BlueChannel) != 0)
00633         {
00634           mean+=p->blue;
00635           sum_squares+=(double) p->blue*p->blue;
00636           sum_cubes+=(double) p->blue*p->blue*p->blue;
00637           sum_fourth_power+=(double) p->blue*p->blue*p->blue*p->blue;
00638           area++;
00639         }
00640       if ((channel & OpacityChannel) != 0)
00641         {
00642           mean+=p->opacity;
00643           sum_squares+=(double) p->opacity*p->opacity;
00644           sum_cubes+=(double) p->opacity*p->opacity*p->opacity;
00645           sum_fourth_power+=(double) p->opacity*p->opacity*p->opacity*
00646             p->opacity;
00647           area++;
00648         }
00649       if (((channel & IndexChannel) != 0) &&
00650           (image->colorspace == CMYKColorspace))
00651         {
00652           mean+=indexes[x];
00653           sum_squares+=(double) indexes[x]*indexes[x];
00654           sum_cubes+=(double) indexes[x]*indexes[x]*indexes[x];
00655           sum_fourth_power+=(double) indexes[x]*indexes[x]*indexes[x]*
00656             indexes[x];
00657           area++;
00658         }
00659       p++;
00660     }
00661   }
00662   if (y < (long) image->rows)
00663     return(MagickFalse);
00664   if (area != 0.0)
00665     {
00666       mean/=area;
00667       sum_squares/=area;
00668       sum_cubes/=area;
00669       sum_fourth_power/=area;
00670     }
00671   standard_deviation=sqrt(sum_squares-(mean*mean));
00672   if (standard_deviation != 0.0)
00673     {
00674       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
00675         3.0*mean*mean*mean*mean;
00676       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
00677         standard_deviation;
00678       *kurtosis-=3.0;
00679       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
00680       *skewness/=standard_deviation*standard_deviation*standard_deviation;
00681     }
00682   return(y == (long) image->rows ? MagickTrue : MagickFalse);
00683 }
00684 
00685 /*
00686 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00687 %                                                                             %
00688 %                                                                             %
00689 %                                                                             %
00690 %   G e t I m a g e C h a n n e l R a n g e                                   %
00691 %                                                                             %
00692 %                                                                             %
00693 %                                                                             %
00694 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00695 %
00696 %  GetImageChannelRange() returns the range of one or more image channels.
00697 %
00698 %  The format of the GetImageChannelRange method is:
00699 %
00700 %      MagickBooleanType GetImageChannelRange(const Image *image,
00701 %        const ChannelType channel,double *minima,double *maxima,
00702 %        ExceptionInfo *exception)
00703 %
00704 %  A description of each parameter follows:
00705 %
00706 %    o image: the image.
00707 %
00708 %    o channel: the channel.
00709 %
00710 %    o minima: the minimum value in the channel.
00711 %
00712 %    o maxima: the maximum value in the channel.
00713 %
00714 %    o exception: return any errors or warnings in this structure.
00715 %
00716 */
00717 
00718 MagickExport MagickBooleanType GetImageRange(const Image *image,
00719   double *minima,double *maxima,ExceptionInfo *exception)
00720 {
00721   return(GetImageChannelRange(image,AllChannels,minima,maxima,exception));
00722 }
00723 
00724 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
00725   const ChannelType channel,double *minima,double *maxima,
00726   ExceptionInfo *exception)
00727 {
00728   long
00729     y;
00730 
00731   MagickPixelPacket
00732     pixel;
00733 
00734   assert(image != (Image *) NULL);
00735   assert(image->signature == MagickSignature);
00736   if (image->debug != MagickFalse)
00737     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00738   *maxima=(-1.0E-37);
00739   *minima=1.0E+37;
00740   GetMagickPixelPacket(image,&pixel);
00741   for (y=0; y < (long) image->rows; y++)
00742   {
00743     register const IndexPacket
00744       *__restrict indexes;
00745 
00746     register const PixelPacket
00747       *__restrict p;
00748 
00749     register long
00750       x;
00751 
00752     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
00753     if (p == (const PixelPacket *) NULL)
00754       break;
00755     indexes=GetVirtualIndexQueue(image);
00756     for (x=0; x < (long) image->columns; x++)
00757     {
00758       SetMagickPixelPacket(image,p,indexes+x,&pixel);
00759       if ((channel & RedChannel) != 0)
00760         {
00761           if (pixel.red < *minima)
00762             *minima=(double) pixel.red;
00763           if (pixel.red > *maxima)
00764             *maxima=(double) pixel.red;
00765         }
00766       if ((channel & GreenChannel) != 0)
00767         {
00768           if (pixel.green < *minima)
00769             *minima=(double) pixel.green;
00770           if (pixel.green > *maxima)
00771             *maxima=(double) pixel.green;
00772         }
00773       if ((channel & BlueChannel) != 0)
00774         {
00775           if (pixel.blue < *minima)
00776             *minima=(double) pixel.blue;
00777           if (pixel.blue > *maxima)
00778             *maxima=(double) pixel.blue;
00779         }
00780       if ((channel & OpacityChannel) != 0)
00781         {
00782           if (pixel.opacity < *minima)
00783             *minima=(double) pixel.opacity;
00784           if (pixel.opacity > *maxima)
00785             *maxima=(double) pixel.opacity;
00786         }
00787       if (((channel & IndexChannel) != 0) &&
00788           (image->colorspace == CMYKColorspace))
00789         {
00790           if ((double) indexes[x] < *minima)
00791             *minima=(double) indexes[x];
00792           if ((double) indexes[x] > *maxima)
00793             *maxima=(double) indexes[x];
00794         }
00795       p++;
00796     }
00797   }
00798   return(y == (long) image->rows ? MagickTrue : MagickFalse);
00799 }
00800 
00801 /*
00802 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00803 %                                                                             %
00804 %                                                                             %
00805 %                                                                             %
00806 %   G e t I m a g e C h a n n e l S t a t i s t i c s                         %
00807 %                                                                             %
00808 %                                                                             %
00809 %                                                                             %
00810 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00811 %
00812 %  GetImageChannelStatistics() returns statistics for each channel in the
00813 %  image.  The statistics include the channel depth, its minima, maxima, mean,
00814 %  standard deviation, kurtosis and skewness.  You can access the red channel
00815 %  mean, for example, like this:
00816 %
00817 %      channel_statistics=GetImageChannelStatistics(image,excepton);
00818 %      red_mean=channel_statistics[RedChannel].mean;
00819 %
00820 %  Use MagickRelinquishMemory() to free the statistics buffer.
00821 %
00822 %  The format of the GetImageChannelStatistics method is:
00823 %
00824 %      ChannelStatistics *GetImageChannelStatistics(const Image *image,
00825 %        ExceptionInfo *exception)
00826 %
00827 %  A description of each parameter follows:
00828 %
00829 %    o image: the image.
00830 %
00831 %    o exception: return any errors or warnings in this structure.
00832 %
00833 */
00834 
00835 static inline double MagickMax(const double x,const double y)
00836 {
00837   if (x > y)
00838     return(x);
00839   return(y);
00840 }
00841 
00842 static inline double MagickMin(const double x,const double y)
00843 {
00844   if (x < y)
00845     return(x);
00846   return(y);
00847 }
00848 
00849 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
00850   ExceptionInfo *exception)
00851 {
00852   ChannelStatistics
00853     *channel_statistics;
00854 
00855   double
00856     area,
00857     sum_squares,
00858     sum_cubes;
00859 
00860   long
00861     y;
00862 
00863   MagickStatusType
00864     status;
00865 
00866   QuantumAny
00867     range;
00868 
00869   register long
00870     i;
00871 
00872   size_t
00873     length;
00874 
00875   unsigned long
00876     channels,
00877     depth;
00878 
00879   assert(image != (Image *) NULL);
00880   assert(image->signature == MagickSignature);
00881   if (image->debug != MagickFalse)
00882     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00883   length=AllChannels+1UL;
00884   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
00885     sizeof(*channel_statistics));
00886   if (channel_statistics == (ChannelStatistics *) NULL)
00887     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
00888   (void) ResetMagickMemory(channel_statistics,0,length*
00889     sizeof(*channel_statistics));
00890   for (i=0; i <= AllChannels; i++)
00891   {
00892     channel_statistics[i].depth=1;
00893     channel_statistics[i].maxima=(-1.0E-37);
00894     channel_statistics[i].minima=1.0E+37;
00895     channel_statistics[i].mean=0.0;
00896     channel_statistics[i].standard_deviation=0.0;
00897     channel_statistics[i].kurtosis=0.0;
00898     channel_statistics[i].skewness=0.0;
00899   }
00900   for (y=0; y < (long) image->rows; y++)
00901   {
00902     register const IndexPacket
00903       *__restrict indexes;
00904 
00905     register const PixelPacket
00906       *__restrict p;
00907 
00908     register long
00909       x;
00910 
00911     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
00912     if (p == (const PixelPacket *) NULL)
00913       break;
00914     indexes=GetVirtualIndexQueue(image);
00915     for (x=0; x < (long) image->columns; )
00916     {
00917       if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
00918         {
00919           depth=channel_statistics[RedChannel].depth;
00920           range=GetQuantumRange(depth);
00921           status=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),
00922             range) ? MagickTrue : MagickFalse;
00923           if (status != MagickFalse)
00924             {
00925               channel_statistics[RedChannel].depth++;
00926               continue;
00927             }
00928         }
00929       if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
00930         {
00931           depth=channel_statistics[GreenChannel].depth;
00932           range=GetQuantumRange(depth);
00933           status=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
00934             range),range) ? MagickTrue : MagickFalse;
00935           if (status != MagickFalse)
00936             {
00937               channel_statistics[GreenChannel].depth++;
00938               continue;
00939             }
00940         }
00941       if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
00942         {
00943           depth=channel_statistics[BlueChannel].depth;
00944           range=GetQuantumRange(depth);
00945           status=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
00946             range),range) ? MagickTrue : MagickFalse;
00947           if (status != MagickFalse)
00948             {
00949               channel_statistics[BlueChannel].depth++;
00950               continue;
00951             }
00952         }
00953       if (image->matte != MagickFalse)
00954         {
00955           if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
00956             {
00957               depth=channel_statistics[OpacityChannel].depth;
00958               range=GetQuantumRange(depth);
00959               status=p->opacity != ScaleAnyToQuantum(ScaleQuantumToAny(
00960                 p->opacity,range),range) ? MagickTrue : MagickFalse;
00961               if (status != MagickFalse)
00962                 {
00963                   channel_statistics[OpacityChannel].depth++;
00964                   continue;
00965                 }
00966             }
00967           }
00968       if (image->colorspace == CMYKColorspace)
00969         {
00970           if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
00971             {
00972               depth=channel_statistics[BlackChannel].depth;
00973               range=GetQuantumRange(depth);
00974               status=indexes[x] != ScaleAnyToQuantum(ScaleQuantumToAny(
00975                 indexes[x],range),range) ? MagickTrue : MagickFalse;
00976               if (status != MagickFalse)
00977                 {
00978                   channel_statistics[BlackChannel].depth++;
00979                   continue;
00980                 }
00981             }
00982         }
00983       if ((double) p->red < channel_statistics[RedChannel].minima)
00984         channel_statistics[RedChannel].minima=(double) p->red;
00985       if ((double) p->red > channel_statistics[RedChannel].maxima)
00986         channel_statistics[RedChannel].maxima=(double) p->red;
00987       channel_statistics[RedChannel].mean+=p->red;
00988       channel_statistics[RedChannel].standard_deviation+=(double) p->red*p->red;
00989       channel_statistics[RedChannel].kurtosis+=(double) p->red*p->red*
00990         p->red*p->red;
00991       channel_statistics[RedChannel].skewness+=(double) p->red*p->red*p->red;
00992       if ((double) p->green < channel_statistics[GreenChannel].minima)
00993         channel_statistics[GreenChannel].minima=(double) p->green;
00994       if ((double) p->green > channel_statistics[GreenChannel].maxima)
00995         channel_statistics[GreenChannel].maxima=(double) p->green;
00996       channel_statistics[GreenChannel].mean+=p->green;
00997       channel_statistics[GreenChannel].standard_deviation+=(double) p->green*
00998         p->green;
00999       channel_statistics[GreenChannel].kurtosis+=(double) p->green*p->green*
01000         p->green*p->green;
01001       channel_statistics[GreenChannel].skewness+=(double) p->green*p->green*
01002         p->green;
01003       if ((double) p->blue < channel_statistics[BlueChannel].minima)
01004         channel_statistics[BlueChannel].minima=(double) p->blue;
01005       if ((double) p->blue > channel_statistics[BlueChannel].maxima)
01006         channel_statistics[BlueChannel].maxima=(double) p->blue;
01007       channel_statistics[BlueChannel].mean+=p->blue;
01008       channel_statistics[BlueChannel].standard_deviation+=(double) p->blue*
01009         p->blue;
01010       channel_statistics[BlueChannel].kurtosis+=(double) p->blue*p->blue*
01011         p->blue*p->blue;
01012       channel_statistics[BlueChannel].skewness+=(double) p->blue*p->blue*
01013         p->blue;
01014       if (image->matte != MagickFalse)
01015         {
01016           if ((double) p->opacity < channel_statistics[OpacityChannel].minima)
01017             channel_statistics[OpacityChannel].minima=(double) p->opacity;
01018           if ((double) p->opacity > channel_statistics[OpacityChannel].maxima)
01019             channel_statistics[OpacityChannel].maxima=(double) p->opacity;
01020           channel_statistics[OpacityChannel].mean+=p->opacity;
01021           channel_statistics[OpacityChannel].standard_deviation+=(double)
01022             p->opacity*p->opacity;
01023           channel_statistics[OpacityChannel].kurtosis+=(double) p->opacity*
01024             p->opacity*p->opacity*p->opacity;
01025           channel_statistics[OpacityChannel].skewness+=(double) p->opacity*
01026             p->opacity*p->opacity;
01027         }
01028       if (image->colorspace == CMYKColorspace)
01029         {
01030           if ((double) indexes[x] < channel_statistics[BlackChannel].minima)
01031             channel_statistics[BlackChannel].minima=(double) indexes[x];
01032           if ((double) indexes[x] > channel_statistics[BlackChannel].maxima)
01033             channel_statistics[BlackChannel].maxima=(double) indexes[x];
01034           channel_statistics[BlackChannel].mean+=indexes[x];
01035           channel_statistics[BlackChannel].standard_deviation+=(double)
01036             indexes[x]*indexes[x];
01037           channel_statistics[BlackChannel].kurtosis+=(double) indexes[x]*
01038             indexes[x]*indexes[x]*indexes[x];
01039           channel_statistics[BlackChannel].skewness+=(double) indexes[x]*
01040             indexes[x]*indexes[x];
01041         }
01042       x++;
01043       p++;
01044     }
01045   }
01046   area=(double) image->columns*image->rows;
01047   for (i=0; i < AllChannels; i++)
01048   {
01049     channel_statistics[i].mean/=area;
01050     channel_statistics[i].standard_deviation/=area;
01051     channel_statistics[i].kurtosis/=area;
01052     channel_statistics[i].skewness/=area;
01053   }
01054   for (i=0; i < AllChannels; i++)
01055   {
01056     channel_statistics[AllChannels].depth=(unsigned long) MagickMax((double)
01057       channel_statistics[AllChannels].depth,(double)
01058       channel_statistics[i].depth);
01059     channel_statistics[AllChannels].minima=MagickMin(
01060       channel_statistics[AllChannels].minima,channel_statistics[i].minima);
01061     channel_statistics[AllChannels].maxima=MagickMax(
01062       channel_statistics[AllChannels].maxima,channel_statistics[i].maxima);
01063     channel_statistics[AllChannels].mean+=channel_statistics[i].mean;
01064     channel_statistics[AllChannels].standard_deviation+=
01065       channel_statistics[i].standard_deviation;
01066     channel_statistics[AllChannels].kurtosis+=channel_statistics[i].kurtosis;
01067     channel_statistics[AllChannels].skewness+=channel_statistics[i].skewness;
01068   }
01069   channels=4;
01070   if (image->colorspace == CMYKColorspace)
01071     channels++;
01072   channel_statistics[AllChannels].mean/=channels;
01073   channel_statistics[AllChannels].standard_deviation/=channels;
01074   channel_statistics[AllChannels].kurtosis/=channels;
01075   channel_statistics[AllChannels].skewness/=channels;
01076   for (i=0; i <= AllChannels; i++)
01077   {
01078     sum_squares=0.0;
01079     sum_squares=channel_statistics[i].standard_deviation;
01080     sum_cubes=0.0;
01081     sum_cubes=channel_statistics[i].skewness;
01082     channel_statistics[i].standard_deviation=sqrt(
01083       channel_statistics[i].standard_deviation-
01084        (channel_statistics[i].mean*channel_statistics[i].mean));
01085     if (channel_statistics[i].standard_deviation == 0.0)
01086       {
01087         channel_statistics[i].kurtosis=0.0;
01088         channel_statistics[i].skewness=0.0;
01089       }
01090     else
01091       {
01092         channel_statistics[i].skewness=(channel_statistics[i].skewness-
01093           3.0*channel_statistics[i].mean*sum_squares+
01094           2.0*channel_statistics[i].mean*channel_statistics[i].mean*
01095           channel_statistics[i].mean)/
01096           (channel_statistics[i].standard_deviation*
01097            channel_statistics[i].standard_deviation*
01098            channel_statistics[i].standard_deviation);
01099         channel_statistics[i].kurtosis=(channel_statistics[i].kurtosis-
01100           4.0*channel_statistics[i].mean*sum_cubes+
01101           6.0*channel_statistics[i].mean*channel_statistics[i].mean*sum_squares-
01102           3.0*channel_statistics[i].mean*channel_statistics[i].mean*
01103           1.0*channel_statistics[i].mean*channel_statistics[i].mean)/
01104           (channel_statistics[i].standard_deviation*
01105            channel_statistics[i].standard_deviation*
01106            channel_statistics[i].standard_deviation*
01107            channel_statistics[i].standard_deviation)-3.0;
01108       }
01109   }
01110   return(channel_statistics);
01111 }

Generated on 19 Nov 2009 for MagickCore by  doxygen 1.6.1