granicus.if.org Git - imagemagick/blob - MagickCore/feature.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %               FFFFF  EEEEE   AAA   TTTTT  U   U  RRRR   EEEEE               %
   7 %               F      E      A   A    T    U   U  R   R  E                   %
   8 %               FFF    EEE    AAAAA    T    U   U  RRRR   EEE                 %
   9 %               F      E      A   A    T    U   U  R R    E                   %
  10 %               F      EEEEE  A   A    T     UUU   R  R   EEEEE               %
  11 %                                                                             %
  12 %                                                                             %
  13 %                      MagickCore Image Feature Methods                       %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                                   Cristy                                    %
  17 %                                 July 1992                                   %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "MagickCore/studio.h"
  44 #include "MagickCore/animate.h"
  45 #include "MagickCore/artifact.h"
  46 #include "MagickCore/blob.h"
  47 #include "MagickCore/blob-private.h"
  48 #include "MagickCore/cache.h"
  49 #include "MagickCore/cache-private.h"
  50 #include "MagickCore/cache-view.h"
  51 #include "MagickCore/client.h"
  52 #include "MagickCore/color.h"
  53 #include "MagickCore/color-private.h"
  54 #include "MagickCore/colorspace.h"
  55 #include "MagickCore/colorspace-private.h"
  56 #include "MagickCore/composite.h"
  57 #include "MagickCore/composite-private.h"
  58 #include "MagickCore/compress.h"
  59 #include "MagickCore/constitute.h"
  60 #include "MagickCore/display.h"
  61 #include "MagickCore/draw.h"
  62 #include "MagickCore/enhance.h"
  63 #include "MagickCore/exception.h"
  64 #include "MagickCore/exception-private.h"
  65 #include "MagickCore/feature.h"
  66 #include "MagickCore/gem.h"
  67 #include "MagickCore/geometry.h"
  68 #include "MagickCore/list.h"
  69 #include "MagickCore/image-private.h"
  70 #include "MagickCore/magic.h"
  71 #include "MagickCore/magick.h"
  72 #include "MagickCore/matrix.h"
  73 #include "MagickCore/memory_.h"
  74 #include "MagickCore/module.h"
  75 #include "MagickCore/monitor.h"
  76 #include "MagickCore/monitor-private.h"
  77 #include "MagickCore/morphology-private.h"
  78 #include "MagickCore/option.h"
  79 #include "MagickCore/paint.h"
  80 #include "MagickCore/pixel-accessor.h"
  81 #include "MagickCore/profile.h"
  82 #include "MagickCore/property.h"
  83 #include "MagickCore/quantize.h"
  84 #include "MagickCore/quantum-private.h"
  85 #include "MagickCore/random_.h"
  86 #include "MagickCore/resource_.h"
  87 #include "MagickCore/segment.h"
  88 #include "MagickCore/semaphore.h"
  89 #include "MagickCore/signature-private.h"
  90 #include "MagickCore/string_.h"
  91 #include "MagickCore/thread-private.h"
  92 #include "MagickCore/timer.h"
  93 #include "MagickCore/utility.h"
  94 #include "MagickCore/version.h"
  95 \f
  96 /*
  97 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  98 %                                                                             %
  99 %                                                                             %
 100 %                                                                             %
 101 %     C a n n y E d g e I m a g e                                             %
 102 %                                                                             %
 103 %                                                                             %
 104 %                                                                             %
 105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 106 %
 107 %  CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
 108 %  edges in images.
 109 %
 110 %  The format of the EdgeImage method is:
 111 %
 112 %      Image *CannyEdgeImage(const Image *image,const double radius,
 113 %        const double sigma,const double lower_percent,
 114 %        const double upper_percent,ExceptionInfo *exception)
 115 %
 116 %  A description of each parameter follows:
 117 %
 118 %    o image: the image.
 119 %
 120 %    o radius: the radius of the gaussian smoothing filter.
 121 %
 122 %    o sigma: the sigma of the gaussian smoothing filter.
 123 %
 124 %    o lower_precent: percentage of edge pixels in the lower threshold.
 125 %
 126 %    o upper_percent: percentage of edge pixels in the upper threshold.
 127 %
 128 %    o exception: return any errors or warnings in this structure.
 129 %
 130 */
 131
 132 typedef struct _CannyInfo
 133 {
 134   double
 135     magnitude,
 136     intensity;
 137
 138   int
 139     orientation;
 140
 141   ssize_t
 142     x,
 143     y;
 144 } CannyInfo;
 145
 146 static inline MagickBooleanType IsAuthenticPixel(const Image *image,
 147   const ssize_t x,const ssize_t y)
 148 {
 149   if ((x < 0) || (x >= (ssize_t) image->columns))
 150     return(MagickFalse);
 151   if ((y < 0) || (y >= (ssize_t) image->rows))
 152     return(MagickFalse);
 153   return(MagickTrue);
 154 }
 155
 156 static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view,
 157   MatrixInfo *canny_cache,const ssize_t x,const ssize_t y,
 158   const double lower_threshold,ExceptionInfo *exception)
 159 {
 160   CannyInfo
 161     edge,
 162     pixel;
 163
 164   MagickBooleanType
 165     status;
 166
 167   register Quantum
 168     *q;
 169
 170   register ssize_t
 171     i;
 172
 173   q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
 174   if (q == (Quantum *) NULL)
 175     return(MagickFalse);
 176   *q=QuantumRange;
 177   status=SyncCacheViewAuthenticPixels(edge_view,exception);
 178   if (status == MagickFalse)
 179     return(MagickFalse);;
 180   if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
 181     return(MagickFalse);
 182   edge.x=x;
 183   edge.y=y;
 184   if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
 185     return(MagickFalse);
 186   for (i=1; i != 0; )
 187   {
 188     ssize_t
 189       v;
 190
 191     i--;
 192     status=GetMatrixElement(canny_cache,i,0,&edge);
 193     if (status == MagickFalse)
 194       return(MagickFalse);
 195     for (v=(-1); v <= 1; v++)
 196     {
 197       ssize_t
 198         u;
 199
 200       for (u=(-1); u <= 1; u++)
 201       {
 202         if ((u == 0) && (v == 0))
 203           continue;
 204         if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
 205           continue;
 206         /*
 207           Not an edge if gradient value is below the lower threshold.
 208         */
 209         q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
 210           exception);
 211         if (q == (Quantum *) NULL)
 212           return(MagickFalse);
 213         status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
 214         if (status == MagickFalse)
 215           return(MagickFalse);
 216         if ((GetPixelIntensity(edge_image,q) == 0.0) &&
 217             (pixel.intensity >= lower_threshold))
 218           {
 219             *q=QuantumRange;
 220             status=SyncCacheViewAuthenticPixels(edge_view,exception);
 221             if (status == MagickFalse)
 222               return(MagickFalse);
 223             edge.x+=u;
 224             edge.y+=v;
 225             status=SetMatrixElement(canny_cache,i,0,&edge);
 226             if (status == MagickFalse)
 227               return(MagickFalse);
 228             i++;
 229           }
 230       }
 231     }
 232   }
 233   return(MagickTrue);
 234 }
 235
 236 MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
 237   const double sigma,const double lower_percent,const double upper_percent,
 238   ExceptionInfo *exception)
 239 {
 240   CacheView
 241     *edge_view;
 242
 243   CannyInfo
 244     pixel;
 245
 246   char
 247     geometry[MaxTextExtent];
 248
 249   double
 250     lower_threshold,
 251     max,
 252     min,
 253     upper_threshold;
 254
 255   Image
 256     *edge_image;
 257
 258   KernelInfo
 259     *kernel_info;
 260
 261   MagickBooleanType
 262     status;
 263
 264   MatrixInfo
 265     *canny_cache;
 266
 267   ssize_t
 268     y;
 269
 270   assert(image != (const Image *) NULL);
 271   assert(image->signature == MagickSignature);
 272   if (image->debug != MagickFalse)
 273     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 274   assert(exception != (ExceptionInfo *) NULL);
 275   assert(exception->signature == MagickSignature);
 276   /*
 277     Filter out noise.
 278   */
 279   (void) FormatLocaleString(geometry,MaxTextExtent,
 280     "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
 281   kernel_info=AcquireKernelInfo(geometry);
 282   if (kernel_info == (KernelInfo *) NULL)
 283     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 284   edge_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
 285     UndefinedCompositeOp,0.0,exception);
 286   kernel_info=DestroyKernelInfo(kernel_info);
 287   if (edge_image == (Image *) NULL)
 288     return((Image *) NULL);
 289   if (SetImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
 290     {
 291       edge_image=DestroyImage(edge_image);
 292       return((Image *) NULL);
 293     }
 294   /*
 295     Find the intensity gradient of the image.
 296   */
 297   canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
 298     sizeof(CannyInfo),exception);
 299   if (canny_cache == (MatrixInfo *) NULL)
 300     {
 301       edge_image=DestroyImage(edge_image);
 302       return((Image *) NULL);
 303     }
 304   status=MagickTrue;
 305   edge_view=AcquireVirtualCacheView(edge_image,exception);
 306 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 307   #pragma omp parallel for schedule(static,4) shared(status) \
 308     magick_threads(edge_image,edge_image,edge_image->rows,1)
 309 #endif
 310   for (y=0; y < (ssize_t) edge_image->rows; y++)
 311   {
 312     register const Quantum
 313       *restrict p;
 314
 315     register ssize_t
 316       x;
 317
 318     if (status == MagickFalse)
 319       continue;
 320     p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
 321       exception);
 322     if (p == (const Quantum *) NULL)
 323       {
 324         status=MagickFalse;
 325         continue;
 326       }
 327     for (x=0; x < (ssize_t) edge_image->columns; x++)
 328     {
 329       CannyInfo
 330         pixel;
 331
 332       double
 333         dx,
 334         dy;
 335
 336       register const Quantum
 337         *restrict kernel_pixels;
 338
 339       ssize_t
 340         v;
 341
 342       static double
 343         Gx[2][2] =
 344         {
 345           { -1.0,  +1.0 },
 346           { -1.0,  +1.0 }
 347         },
 348         Gy[2][2] =
 349         {
 350           { +1.0, +1.0 },
 351           { -1.0, -1.0 }
 352         };
 353
 354       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
 355       dx=0.0;
 356       dy=0.0;
 357       kernel_pixels=p;
 358       for (v=0; v < 2; v++)
 359       {
 360         ssize_t
 361           u;
 362
 363         for (u=0; u < 2; u++)
 364         {
 365           double
 366             intensity;
 367
 368           intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
 369           dx+=0.5*Gx[v][u]*intensity;
 370           dy+=0.5*Gy[v][u]*intensity;
 371         }
 372         kernel_pixels+=edge_image->columns+1;
 373       }
 374       pixel.magnitude=hypot(dx,dy);
 375       pixel.orientation=0;
 376       if (fabs(dx) > MagickEpsilon)
 377         {
 378           double
 379             slope;
 380
 381           slope=dy/dx;
 382           if (slope < 0.0)
 383             {
 384               if (slope < -2.41421356237)
 385                 pixel.orientation=0;
 386               else
 387                 if (slope < -0.414213562373)
 388                   pixel.orientation=1;
 389                 else
 390                   pixel.orientation=2;
 391             }
 392           else
 393             {
 394               if (slope > 2.41421356237)
 395                 pixel.orientation=0;
 396               else
 397                 if (slope > 0.414213562373)
 398                   pixel.orientation=3;
 399                 else
 400                   pixel.orientation=2;
 401             }
 402         }
 403       if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
 404         continue;
 405       p+=GetPixelChannels(edge_image);
 406     }
 407   }
 408   edge_view=DestroyCacheView(edge_view);
 409   /*
 410     Non-maxima suppression, remove pixels that are not considered to be part
 411     of an edge.
 412   */
 413   (void) GetMatrixElement(canny_cache,0,0,&pixel);
 414   max=pixel.intensity;
 415   min=pixel.intensity;
 416   edge_view=AcquireAuthenticCacheView(edge_image,exception);
 417 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 418   #pragma omp parallel for schedule(static,4) shared(status) \
 419     magick_threads(edge_image,edge_image,edge_image->rows,1)
 420 #endif
 421   for (y=0; y < (ssize_t) edge_image->rows; y++)
 422   {
 423     register Quantum
 424       *restrict q;
 425
 426     register ssize_t
 427       x;
 428
 429     if (status == MagickFalse)
 430       continue;
 431     q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
 432       exception);
 433     if (q == (Quantum *) NULL)
 434       {
 435         status=MagickFalse;
 436         continue;
 437       }
 438     for (x=0; x < (ssize_t) edge_image->columns; x++)
 439     {
 440       CannyInfo
 441         alpha_pixel,
 442         beta_pixel,
 443         pixel;
 444
 445       (void) GetMatrixElement(canny_cache,x,y,&pixel);
 446       switch (pixel.orientation)
 447       {
 448         case 0:
 449         default:
 450         {
 451           /*
 452             0 degrees, north and south.
 453           */
 454           (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
 455           (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
 456           break;
 457         }
 458         case 1:
 459         {
 460           /*
 461             45 degrees, northwest and southeast.
 462           */
 463           (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
 464           (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
 465           break;
 466         }
 467         case 2:
 468         {
 469           /*
 470             90 degrees, east and west.
 471           */
 472           (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
 473           (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
 474           break;
 475         }
 476         case 3:
 477         {
 478           /*
 479             135 degrees, northeast and southwest.
 480           */
 481           (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
 482           (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
 483           break;
 484         }
 485       }
 486       pixel.intensity=pixel.magnitude;
 487       if ((pixel.magnitude < alpha_pixel.magnitude) ||
 488           (pixel.magnitude < beta_pixel.magnitude))
 489         pixel.intensity=0;
 490       (void) SetMatrixElement(canny_cache,x,y,&pixel);
 491 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 492       #pragma omp critical (MagickCore_CannyEdgeImage)
 493 #endif
 494       {
 495         if (pixel.intensity < min)
 496           min=pixel.intensity;
 497         if (pixel.intensity > max)
 498           max=pixel.intensity;
 499       }
 500       *q=0;
 501       q+=GetPixelChannels(edge_image);
 502     }
 503     if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
 504       status=MagickFalse;
 505   }
 506   edge_view=DestroyCacheView(edge_view);
 507   /*
 508     Estimate hysteresis threshold.
 509   */
 510   lower_threshold=lower_percent*(max-min)+min;
 511   upper_threshold=upper_percent*(max-min)+min;
 512   /*
 513     Hysteresis threshold.
 514   */
 515   edge_view=AcquireAuthenticCacheView(edge_image,exception);
 516   for (y=0; y < (ssize_t) edge_image->rows; y++)
 517   {
 518     register ssize_t
 519       x;
 520
 521     if (status == MagickFalse)
 522       continue;
 523     for (x=0; x < (ssize_t) edge_image->columns; x++)
 524     {
 525       CannyInfo
 526         pixel;
 527
 528       register const Quantum
 529         *restrict p;
 530
 531       /*
 532         Edge if pixel gradient higher than upper threshold.
 533       */
 534       p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
 535       if (p == (const Quantum *) NULL)
 536         continue;
 537       status=GetMatrixElement(canny_cache,x,y,&pixel);
 538       if (status == MagickFalse)
 539         continue;
 540       if ((GetPixelIntensity(edge_image,p) == 0.0) &&
 541           (pixel.intensity >= upper_threshold))
 542         status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
 543           exception);
 544     }
 545   }
 546   edge_view=DestroyCacheView(edge_view);
 547   /*
 548     Free resources.
 549   */
 550   canny_cache=DestroyMatrixInfo(canny_cache);
 551   return(edge_image);
 552 }
 553 \f
 554 /*
 555 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 556 %                                                                             %
 557 %                                                                             %
 558 %                                                                             %
 559 %     H o u g h L i n e s I m a g e                                           %
 560 %                                                                             %
 561 %                                                                             %
 562 %                                                                             %
 563 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 564 %
 565 %  HoughLinesImage() identifies lines in the image.
 566 %
 567 %  The format of the HoughLinesImage method is:
 568 %
 569 %      Image *HoughLinesImage(const Image *image,const size_t width,
 570 %        const size_t height,const size_t threshold,ExceptionInfo *exception)
 571 %
 572 %  A description of each parameter follows:
 573 %
 574 %    o image: the image.
 575 %
 576 %    o width, height: find line pairs as local maxima in this neighborhood.
 577 %
 578 %    o threshold: the line count threshold.
 579 %
 580 %    o exception: return any errors or warnings in this structure.
 581 %
 582 */
 583
 584 static inline double MagickRound(double x)
 585 {
 586   /*
 587     Round the fraction to nearest integer.
 588   */
 589   if ((x-floor(x)) < (ceil(x)-x))
 590     return(floor(x));
 591   return(ceil(x));
 592 }
 593
 594 MagickExport Image *HoughLinesImage(const Image *image,const size_t width,
 595   const size_t height,const size_t threshold,ExceptionInfo *exception)
 596 {
 597   CacheView
 598     *image_view;
 599
 600   char
 601     message[MaxTextExtent],
 602     path[MaxTextExtent];
 603
 604   const char
 605     *artifact;
 606
 607   double
 608     hough_height;
 609
 610   Image
 611     *lines_image = NULL;
 612
 613   ImageInfo
 614     *image_info;
 615
 616   int
 617     file;
 618
 619   MagickBooleanType
 620     status;
 621
 622   MatrixInfo
 623     *accumulator;
 624
 625   PointInfo
 626     center;
 627
 628   register ssize_t
 629     y;
 630
 631   size_t
 632     accumulator_height,
 633     accumulator_width,
 634     line_count;
 635
 636   /*
 637     Create the accumulator.
 638   */
 639   assert(image != (const Image *) NULL);
 640   assert(image->signature == MagickSignature);
 641   if (image->debug != MagickFalse)
 642     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 643   assert(exception != (ExceptionInfo *) NULL);
 644   assert(exception->signature == MagickSignature);
 645   accumulator_width=180;
 646   hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
 647     image->rows : image->columns))/2.0);
 648   accumulator_height=(size_t) (2.0*hough_height);
 649   accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
 650     sizeof(double),exception);
 651   if (accumulator == (MatrixInfo *) NULL)
 652     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 653   if (NullMatrix(accumulator) == MagickFalse)
 654     {
 655       accumulator=DestroyMatrixInfo(accumulator);
 656       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 657     }
 658   /*
 659     Populate the accumulator.
 660   */
 661   status=MagickTrue;
 662   center.x=(double) image->columns/2.0;
 663   center.y=(double) image->rows/2.0;
 664   image_view=AcquireVirtualCacheView(image,exception);
 665   for (y=0; y < (ssize_t) image->rows; y++)
 666   {
 667     register const Quantum
 668       *restrict p;
 669
 670     register ssize_t
 671       x;
 672
 673     if (status == MagickFalse)
 674       continue;
 675     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 676     if (p == (Quantum *) NULL)
 677       {
 678         status=MagickFalse;
 679         continue;
 680       }
 681     for (x=0; x < (ssize_t) image->columns; x++)
 682     {
 683       if (GetPixelIntensity(image,p) > (QuantumRange/2))
 684         {
 685           register ssize_t
 686             i;
 687
 688           for (i=0; i < 180; i++)
 689           {
 690             double
 691               count,
 692               radius;
 693
 694             radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+
 695               (((double) y-center.y)*sin(DegreesToRadians((double) i)));
 696             (void) GetMatrixElement(accumulator,i,(ssize_t)
 697               MagickRound(radius+hough_height),&count);
 698             count++;
 699             (void) SetMatrixElement(accumulator,i,(ssize_t)
 700               MagickRound(radius+hough_height),&count);
 701           }
 702         }
 703       p+=GetPixelChannels(image);
 704     }
 705   }
 706   image_view=DestroyCacheView(image_view);
 707   if (status == MagickFalse)
 708     {
 709       accumulator=DestroyMatrixInfo(accumulator);
 710       return((Image *) NULL);
 711     }
 712   /*
 713     Generate line segments from accumulator.
 714   */
 715   file=AcquireUniqueFileResource(path);
 716   if (file == -1)
 717     {
 718       accumulator=DestroyMatrixInfo(accumulator);
 719       return((Image *) NULL);
 720     }
 721   (void) FormatLocaleString(message,MaxTextExtent,
 722     "# Hough line transform: %.20gx%.20g%+.20g\n",(double) width,
 723     (double) height,(double) threshold);
 724   (void) write(file,message,strlen(message));
 725   (void) FormatLocaleString(message,MaxTextExtent,"viewbox 0 0 %.20g %.20g\n",
 726     (double) image->columns,(double) image->rows);
 727   (void) write(file,message,strlen(message));
 728   line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
 729   if (threshold != 0)
 730     line_count=threshold;
 731   for (y=0; y < (ssize_t) accumulator_height; y++)
 732   {
 733     register ssize_t
 734       x;
 735
 736     for (x=0; x < (ssize_t) accumulator_width; x++)
 737     {
 738       double
 739         count;
 740
 741       (void) GetMatrixElement(accumulator,x,y,&count);
 742       if (count >= (double) line_count)
 743         {
 744           double
 745             maxima;
 746
 747           SegmentInfo
 748             line;
 749
 750           ssize_t
 751             v;
 752
 753           /*
 754             Is point a local maxima?
 755           */
 756           maxima=count;
 757           for (v=(-((ssize_t) height/2)); v < (((ssize_t) height/2)); v++)
 758           {
 759             ssize_t
 760               u;
 761
 762             for (u=(-((ssize_t) width/2)); u < (((ssize_t) width/2)); u++)
 763             {
 764               if ((u != 0) || (v !=0))
 765                 {
 766                   (void) GetMatrixElement(accumulator,x+u,y+v,&count);
 767                   if (count > maxima)
 768                     {
 769                       maxima=count;
 770                       break;
 771                     }
 772                 }
 773             }
 774             if (u < (ssize_t) (width/2))
 775               break;
 776           }
 777           (void) GetMatrixElement(accumulator,x,y,&count);
 778           if (maxima > count)
 779             continue;
 780           if ((x >= 45) && (x <= 135))
 781             {
 782               /*
 783                 y = (r-x cos(t))/sin(t)
 784               */
 785               line.x1=0.0;
 786               line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
 787                 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
 788                 sin(DegreesToRadians((double) x))+(image->rows/2.0);
 789               line.x2=(double) image->columns;
 790               line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
 791                 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
 792                 sin(DegreesToRadians((double) x))+(image->rows/2.0);
 793             }
 794           else
 795             {
 796               /*
 797                 x = (r-y cos(t))/sin(t)
 798               */
 799               line.y1=0.0;
 800               line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
 801                 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
 802                 cos(DegreesToRadians((double) x))+(image->columns/2.0);
 803               line.y2=(double) image->rows;
 804               line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
 805                 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
 806                 cos(DegreesToRadians((double) x))+(image->columns/2.0);
 807             }
 808           (void) FormatLocaleString(message,MaxTextExtent,
 809             "line %g,%g %g,%g  # %g\n",line.x1,line.y1,line.x2,line.y2,maxima);
 810           (void) write(file,message,strlen(message));
 811         }
 812     }
 813   }
 814   (void) close(file);
 815   /*
 816     Render lines to image canvas.
 817   */
 818   image_info=AcquireImageInfo();
 819   image_info->background_color=image->background_color;
 820   (void) FormatLocaleString(image_info->filename,MaxTextExtent,"mvg:%s",path);
 821   artifact=GetImageArtifact(image,"background");
 822   if (artifact != (const char *) NULL)
 823     (void) SetImageOption(image_info,"background",artifact);
 824   artifact=GetImageArtifact(image,"fill");
 825   if (artifact != (const char *) NULL)
 826     (void) SetImageOption(image_info,"fill",artifact);
 827   artifact=GetImageArtifact(image,"stroke");
 828   if (artifact != (const char *) NULL)
 829     (void) SetImageOption(image_info,"stroke",artifact);
 830   artifact=GetImageArtifact(image,"strokewidth");
 831   if (artifact != (const char *) NULL)
 832     (void) SetImageOption(image_info,"strokewidth",artifact);
 833   lines_image=ReadImage(image_info,exception);
 834   artifact=GetImageArtifact(image,"hough-lines:accumulator");
 835   if ((lines_image != (Image *) NULL) &&
 836       (IsStringTrue(artifact) != MagickFalse))
 837     {
 838       Image
 839         *accumulator_image;
 840
 841       accumulator_image=MatrixToImage(accumulator,exception);
 842       if (accumulator_image != (Image *) NULL)
 843         AppendImageToList(&lines_image,accumulator_image);
 844     }
 845   /*
 846     Free resources.
 847   */
 848   accumulator=DestroyMatrixInfo(accumulator);
 849   image_info=DestroyImageInfo(image_info);
 850   (void) RelinquishUniqueFileResource(path);
 851   return(GetFirstImageInList(lines_image));
 852 }
 853 \f
 854 /*
 855 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 856 %                                                                             %
 857 %                                                                             %
 858 %                                                                             %
 859 %   G e t I m a g e F e a t u r e s                                           %
 860 %                                                                             %
 861 %                                                                             %
 862 %                                                                             %
 863 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 864 %
 865 %  GetImageFeatures() returns features for each channel in the image in
 866 %  each of four directions (horizontal, vertical, left and right diagonals)
 867 %  for the specified distance.  The features include the angular second
 868 %  moment, contrast, correlation, sum of squares: variance, inverse difference
 869 %  moment, sum average, sum varience, sum entropy, entropy, difference variance,%  difference entropy, information measures of correlation 1, information
 870 %  measures of correlation 2, and maximum correlation coefficient.  You can
 871 %  access the red channel contrast, for example, like this:
 872 %
 873 %      channel_features=GetImageFeatures(image,1,exception);
 874 %      contrast=channel_features[RedPixelChannel].contrast[0];
 875 %
 876 %  Use MagickRelinquishMemory() to free the features buffer.
 877 %
 878 %  The format of the GetImageFeatures method is:
 879 %
 880 %      ChannelFeatures *GetImageFeatures(const Image *image,
 881 %        const size_t distance,ExceptionInfo *exception)
 882 %
 883 %  A description of each parameter follows:
 884 %
 885 %    o image: the image.
 886 %
 887 %    o distance: the distance.
 888 %
 889 %    o exception: return any errors or warnings in this structure.
 890 %
 891 */
 892
 893 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
 894 {
 895   if (x < 0)
 896     return(-x);
 897   return(x);
 898 }
 899
 900 static inline double MagickLog10(const double x)
 901 {
 902 #define Log10Epsilon  (1.0e-11)
 903
 904  if (fabs(x) < Log10Epsilon)
 905    return(log10(Log10Epsilon));
 906  return(log10(fabs(x)));
 907 }
 908
 909 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
 910   const size_t distance,ExceptionInfo *exception)
 911 {
 912   typedef struct _ChannelStatistics
 913   {
 914     PixelInfo
 915       direction[4];  /* horizontal, vertical, left and right diagonals */
 916   } ChannelStatistics;
 917
 918   CacheView
 919     *image_view;
 920
 921   ChannelFeatures
 922     *channel_features;
 923
 924   ChannelStatistics
 925     **cooccurrence,
 926     correlation,
 927     *density_x,
 928     *density_xy,
 929     *density_y,
 930     entropy_x,
 931     entropy_xy,
 932     entropy_xy1,
 933     entropy_xy2,
 934     entropy_y,
 935     mean,
 936     **Q,
 937     *sum,
 938     sum_squares,
 939     variance;
 940
 941   PixelPacket
 942     gray,
 943     *grays;
 944
 945   MagickBooleanType
 946     status;
 947
 948   register ssize_t
 949     i;
 950
 951   size_t
 952     length;
 953
 954   ssize_t
 955     y;
 956
 957   unsigned int
 958     number_grays;
 959
 960   assert(image != (Image *) NULL);
 961   assert(image->signature == MagickSignature);
 962   if (image->debug != MagickFalse)
 963     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 964   if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
 965     return((ChannelFeatures *) NULL);
 966   length=CompositeChannels+1UL;
 967   channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
 968     sizeof(*channel_features));
 969   if (channel_features == (ChannelFeatures *) NULL)
 970     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 971   (void) ResetMagickMemory(channel_features,0,length*
 972     sizeof(*channel_features));
 973   /*
 974     Form grays.
 975   */
 976   grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
 977   if (grays == (PixelPacket *) NULL)
 978     {
 979       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 980         channel_features);
 981       (void) ThrowMagickException(exception,GetMagickModule(),
 982         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 983       return(channel_features);
 984     }
 985   for (i=0; i <= (ssize_t) MaxMap; i++)
 986   {
 987     grays[i].red=(~0U);
 988     grays[i].green=(~0U);
 989     grays[i].blue=(~0U);
 990     grays[i].alpha=(~0U);
 991     grays[i].black=(~0U);
 992   }
 993   status=MagickTrue;
 994   image_view=AcquireVirtualCacheView(image,exception);
 995 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 996   #pragma omp parallel for schedule(static,4) shared(status) \
 997     magick_threads(image,image,image->rows,1)
 998 #endif
 999   for (y=0; y < (ssize_t) image->rows; y++)
1000   {
1001     register const Quantum
1002       *restrict p;
1003
1004     register ssize_t
1005       x;
1006
1007     if (status == MagickFalse)
1008       continue;
1009     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1010     if (p == (const Quantum *) NULL)
1011       {
1012         status=MagickFalse;
1013         continue;
1014       }
1015     for (x=0; x < (ssize_t) image->columns; x++)
1016     {
1017       grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
1018         ScaleQuantumToMap(GetPixelRed(image,p));
1019       grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
1020         ScaleQuantumToMap(GetPixelGreen(image,p));
1021       grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
1022         ScaleQuantumToMap(GetPixelBlue(image,p));
1023       if (image->colorspace == CMYKColorspace)
1024         grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
1025           ScaleQuantumToMap(GetPixelBlack(image,p));
1026       if (image->alpha_trait == BlendPixelTrait)
1027         grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
1028           ScaleQuantumToMap(GetPixelAlpha(image,p));
1029       p+=GetPixelChannels(image);
1030     }
1031   }
1032   image_view=DestroyCacheView(image_view);
1033   if (status == MagickFalse)
1034     {
1035       grays=(PixelPacket *) RelinquishMagickMemory(grays);
1036       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1037         channel_features);
1038       return(channel_features);
1039     }
1040   (void) ResetMagickMemory(&gray,0,sizeof(gray));
1041   for (i=0; i <= (ssize_t) MaxMap; i++)
1042   {
1043     if (grays[i].red != ~0U)
1044       grays[gray.red++].red=grays[i].red;
1045     if (grays[i].green != ~0U)
1046       grays[gray.green++].green=grays[i].green;
1047     if (grays[i].blue != ~0U)
1048       grays[gray.blue++].blue=grays[i].blue;
1049     if (image->colorspace == CMYKColorspace)
1050       if (grays[i].black != ~0U)
1051         grays[gray.black++].black=grays[i].black;
1052     if (image->alpha_trait == BlendPixelTrait)
1053       if (grays[i].alpha != ~0U)
1054         grays[gray.alpha++].alpha=grays[i].alpha;
1055   }
1056   /*
1057     Allocate spatial dependence matrix.
1058   */
1059   number_grays=gray.red;
1060   if (gray.green > number_grays)
1061     number_grays=gray.green;
1062   if (gray.blue > number_grays)
1063     number_grays=gray.blue;
1064   if (image->colorspace == CMYKColorspace)
1065     if (gray.black > number_grays)
1066       number_grays=gray.black;
1067   if (image->alpha_trait == BlendPixelTrait)
1068     if (gray.alpha > number_grays)
1069       number_grays=gray.alpha;
1070   cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
1071     sizeof(*cooccurrence));
1072   density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
1073     sizeof(*density_x));
1074   density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
1075     sizeof(*density_xy));
1076   density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
1077     sizeof(*density_y));
1078   Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
1079   sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
1080   if ((cooccurrence == (ChannelStatistics **) NULL) ||
1081       (density_x == (ChannelStatistics *) NULL) ||
1082       (density_xy == (ChannelStatistics *) NULL) ||
1083       (density_y == (ChannelStatistics *) NULL) ||
1084       (Q == (ChannelStatistics **) NULL) ||
1085       (sum == (ChannelStatistics *) NULL))
1086     {
1087       if (Q != (ChannelStatistics **) NULL)
1088         {
1089           for (i=0; i < (ssize_t) number_grays; i++)
1090             Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1091           Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1092         }
1093       if (sum != (ChannelStatistics *) NULL)
1094         sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1095       if (density_y != (ChannelStatistics *) NULL)
1096         density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1097       if (density_xy != (ChannelStatistics *) NULL)
1098         density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1099       if (density_x != (ChannelStatistics *) NULL)
1100         density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1101       if (cooccurrence != (ChannelStatistics **) NULL)
1102         {
1103           for (i=0; i < (ssize_t) number_grays; i++)
1104             cooccurrence[i]=(ChannelStatistics *)
1105               RelinquishMagickMemory(cooccurrence[i]);
1106           cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
1107             cooccurrence);
1108         }
1109       grays=(PixelPacket *) RelinquishMagickMemory(grays);
1110       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1111         channel_features);
1112       (void) ThrowMagickException(exception,GetMagickModule(),
1113         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1114       return(channel_features);
1115     }
1116   (void) ResetMagickMemory(&correlation,0,sizeof(correlation));
1117   (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x));
1118   (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
1119   (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y));
1120   (void) ResetMagickMemory(&mean,0,sizeof(mean));
1121   (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum));
1122   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
1123   (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy));
1124   (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x));
1125   (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy));
1126   (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1));
1127   (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2));
1128   (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y));
1129   (void) ResetMagickMemory(&variance,0,sizeof(variance));
1130   for (i=0; i < (ssize_t) number_grays; i++)
1131   {
1132     cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
1133       sizeof(**cooccurrence));
1134     Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
1135     if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
1136         (Q[i] == (ChannelStatistics *) NULL))
1137       break;
1138     (void) ResetMagickMemory(cooccurrence[i],0,number_grays*
1139       sizeof(**cooccurrence));
1140     (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q));
1141   }
1142   if (i < (ssize_t) number_grays)
1143     {
1144       for (i--; i >= 0; i--)
1145       {
1146         if (Q[i] != (ChannelStatistics *) NULL)
1147           Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1148         if (cooccurrence[i] != (ChannelStatistics *) NULL)
1149           cooccurrence[i]=(ChannelStatistics *)
1150             RelinquishMagickMemory(cooccurrence[i]);
1151       }
1152       Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1153       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1154       sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1155       density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1156       density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1157       density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1158       grays=(PixelPacket *) RelinquishMagickMemory(grays);
1159       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1160         channel_features);
1161       (void) ThrowMagickException(exception,GetMagickModule(),
1162         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1163       return(channel_features);
1164     }
1165   /*
1166     Initialize spatial dependence matrix.
1167   */
1168   status=MagickTrue;
1169   image_view=AcquireVirtualCacheView(image,exception);
1170   for (y=0; y < (ssize_t) image->rows; y++)
1171   {
1172     register const Quantum
1173       *restrict p;
1174
1175     register ssize_t
1176       x;
1177
1178     ssize_t
1179       i,
1180       offset,
1181       u,
1182       v;
1183
1184     if (status == MagickFalse)
1185       continue;
1186     p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+
1187       2*distance,distance+2,exception);
1188     if (p == (const Quantum *) NULL)
1189       {
1190         status=MagickFalse;
1191         continue;
1192       }
1193     p+=distance*GetPixelChannels(image);;
1194     for (x=0; x < (ssize_t) image->columns; x++)
1195     {
1196       for (i=0; i < 4; i++)
1197       {
1198         switch (i)
1199         {
1200           case 0:
1201           default:
1202           {
1203             /*
1204               Horizontal adjacency.
1205             */
1206             offset=(ssize_t) distance;
1207             break;
1208           }
1209           case 1:
1210           {
1211             /*
1212               Vertical adjacency.
1213             */
1214             offset=(ssize_t) (image->columns+2*distance);
1215             break;
1216           }
1217           case 2:
1218           {
1219             /*
1220               Right diagonal adjacency.
1221             */
1222             offset=(ssize_t) ((image->columns+2*distance)-distance);
1223             break;
1224           }
1225           case 3:
1226           {
1227             /*
1228               Left diagonal adjacency.
1229             */
1230             offset=(ssize_t) ((image->columns+2*distance)+distance);
1231             break;
1232           }
1233         }
1234         u=0;
1235         v=0;
1236         while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
1237           u++;
1238         while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
1239           v++;
1240         cooccurrence[u][v].direction[i].red++;
1241         cooccurrence[v][u].direction[i].red++;
1242         u=0;
1243         v=0;
1244         while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
1245           u++;
1246         while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
1247           v++;
1248         cooccurrence[u][v].direction[i].green++;
1249         cooccurrence[v][u].direction[i].green++;
1250         u=0;
1251         v=0;
1252         while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
1253           u++;
1254         while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
1255           v++;
1256         cooccurrence[u][v].direction[i].blue++;
1257         cooccurrence[v][u].direction[i].blue++;
1258         if (image->colorspace == CMYKColorspace)
1259           {
1260             u=0;
1261             v=0;
1262             while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
1263               u++;
1264             while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
1265               v++;
1266             cooccurrence[u][v].direction[i].black++;
1267             cooccurrence[v][u].direction[i].black++;
1268           }
1269         if (image->alpha_trait == BlendPixelTrait)
1270           {
1271             u=0;
1272             v=0;
1273             while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
1274               u++;
1275             while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
1276               v++;
1277             cooccurrence[u][v].direction[i].alpha++;
1278             cooccurrence[v][u].direction[i].alpha++;
1279           }
1280       }
1281       p+=GetPixelChannels(image);
1282     }
1283   }
1284   grays=(PixelPacket *) RelinquishMagickMemory(grays);
1285   image_view=DestroyCacheView(image_view);
1286   if (status == MagickFalse)
1287     {
1288       for (i=0; i < (ssize_t) number_grays; i++)
1289         cooccurrence[i]=(ChannelStatistics *)
1290           RelinquishMagickMemory(cooccurrence[i]);
1291       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1292       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1293         channel_features);
1294       (void) ThrowMagickException(exception,GetMagickModule(),
1295         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1296       return(channel_features);
1297     }
1298   /*
1299     Normalize spatial dependence matrix.
1300   */
1301   for (i=0; i < 4; i++)
1302   {
1303     double
1304       normalize;
1305
1306     register ssize_t
1307       y;
1308
1309     switch (i)
1310     {
1311       case 0:
1312       default:
1313       {
1314         /*
1315           Horizontal adjacency.
1316         */
1317         normalize=2.0*image->rows*(image->columns-distance);
1318         break;
1319       }
1320       case 1:
1321       {
1322         /*
1323           Vertical adjacency.
1324         */
1325         normalize=2.0*(image->rows-distance)*image->columns;
1326         break;
1327       }
1328       case 2:
1329       {
1330         /*
1331           Right diagonal adjacency.
1332         */
1333         normalize=2.0*(image->rows-distance)*(image->columns-distance);
1334         break;
1335       }
1336       case 3:
1337       {
1338         /*
1339           Left diagonal adjacency.
1340         */
1341         normalize=2.0*(image->rows-distance)*(image->columns-distance);
1342         break;
1343       }
1344     }
1345     normalize=PerceptibleReciprocal(normalize);
1346     for (y=0; y < (ssize_t) number_grays; y++)
1347     {
1348       register ssize_t
1349         x;
1350
1351       for (x=0; x < (ssize_t) number_grays; x++)
1352       {
1353         cooccurrence[x][y].direction[i].red*=normalize;
1354         cooccurrence[x][y].direction[i].green*=normalize;
1355         cooccurrence[x][y].direction[i].blue*=normalize;
1356         if (image->colorspace == CMYKColorspace)
1357           cooccurrence[x][y].direction[i].black*=normalize;
1358         if (image->alpha_trait == BlendPixelTrait)
1359           cooccurrence[x][y].direction[i].alpha*=normalize;
1360       }
1361     }
1362   }
1363   /*
1364     Compute texture features.
1365   */
1366 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1367   #pragma omp parallel for schedule(static,4) shared(status) \
1368     magick_threads(image,image,number_grays,1)
1369 #endif
1370   for (i=0; i < 4; i++)
1371   {
1372     register ssize_t
1373       y;
1374
1375     for (y=0; y < (ssize_t) number_grays; y++)
1376     {
1377       register ssize_t
1378         x;
1379
1380       for (x=0; x < (ssize_t) number_grays; x++)
1381       {
1382         /*
1383           Angular second moment:  measure of homogeneity of the image.
1384         */
1385         channel_features[RedPixelChannel].angular_second_moment[i]+=
1386           cooccurrence[x][y].direction[i].red*
1387           cooccurrence[x][y].direction[i].red;
1388         channel_features[GreenPixelChannel].angular_second_moment[i]+=
1389           cooccurrence[x][y].direction[i].green*
1390           cooccurrence[x][y].direction[i].green;
1391         channel_features[BluePixelChannel].angular_second_moment[i]+=
1392           cooccurrence[x][y].direction[i].blue*
1393           cooccurrence[x][y].direction[i].blue;
1394         if (image->colorspace == CMYKColorspace)
1395           channel_features[BlackPixelChannel].angular_second_moment[i]+=
1396             cooccurrence[x][y].direction[i].black*
1397             cooccurrence[x][y].direction[i].black;
1398         if (image->alpha_trait == BlendPixelTrait)
1399           channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1400             cooccurrence[x][y].direction[i].alpha*
1401             cooccurrence[x][y].direction[i].alpha;
1402         /*
1403           Correlation: measure of linear-dependencies in the image.
1404         */
1405         sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1406         sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1407         sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1408         if (image->colorspace == CMYKColorspace)
1409           sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1410         if (image->alpha_trait == BlendPixelTrait)
1411           sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1412         correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1413         correlation.direction[i].green+=x*y*
1414           cooccurrence[x][y].direction[i].green;
1415         correlation.direction[i].blue+=x*y*
1416           cooccurrence[x][y].direction[i].blue;
1417         if (image->colorspace == CMYKColorspace)
1418           correlation.direction[i].black+=x*y*
1419             cooccurrence[x][y].direction[i].black;
1420         if (image->alpha_trait == BlendPixelTrait)
1421           correlation.direction[i].alpha+=x*y*
1422             cooccurrence[x][y].direction[i].alpha;
1423         /*
1424           Inverse Difference Moment.
1425         */
1426         channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1427           cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1428         channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1429           cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1430         channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1431           cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1432         if (image->colorspace == CMYKColorspace)
1433           channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1434             cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1435         if (image->alpha_trait == BlendPixelTrait)
1436           channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1437             cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1438         /*
1439           Sum average.
1440         */
1441         density_xy[y+x+2].direction[i].red+=
1442           cooccurrence[x][y].direction[i].red;
1443         density_xy[y+x+2].direction[i].green+=
1444           cooccurrence[x][y].direction[i].green;
1445         density_xy[y+x+2].direction[i].blue+=
1446           cooccurrence[x][y].direction[i].blue;
1447         if (image->colorspace == CMYKColorspace)
1448           density_xy[y+x+2].direction[i].black+=
1449             cooccurrence[x][y].direction[i].black;
1450         if (image->alpha_trait == BlendPixelTrait)
1451           density_xy[y+x+2].direction[i].alpha+=
1452             cooccurrence[x][y].direction[i].alpha;
1453         /*
1454           Entropy.
1455         */
1456         channel_features[RedPixelChannel].entropy[i]-=
1457           cooccurrence[x][y].direction[i].red*
1458           MagickLog10(cooccurrence[x][y].direction[i].red);
1459         channel_features[GreenPixelChannel].entropy[i]-=
1460           cooccurrence[x][y].direction[i].green*
1461           MagickLog10(cooccurrence[x][y].direction[i].green);
1462         channel_features[BluePixelChannel].entropy[i]-=
1463           cooccurrence[x][y].direction[i].blue*
1464           MagickLog10(cooccurrence[x][y].direction[i].blue);
1465         if (image->colorspace == CMYKColorspace)
1466           channel_features[BlackPixelChannel].entropy[i]-=
1467             cooccurrence[x][y].direction[i].black*
1468             MagickLog10(cooccurrence[x][y].direction[i].black);
1469         if (image->alpha_trait == BlendPixelTrait)
1470           channel_features[AlphaPixelChannel].entropy[i]-=
1471             cooccurrence[x][y].direction[i].alpha*
1472             MagickLog10(cooccurrence[x][y].direction[i].alpha);
1473         /*
1474           Information Measures of Correlation.
1475         */
1476         density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1477         density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1478         density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1479         if (image->alpha_trait == BlendPixelTrait)
1480           density_x[x].direction[i].alpha+=
1481             cooccurrence[x][y].direction[i].alpha;
1482         if (image->colorspace == CMYKColorspace)
1483           density_x[x].direction[i].black+=
1484             cooccurrence[x][y].direction[i].black;
1485         density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1486         density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1487         density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1488         if (image->colorspace == CMYKColorspace)
1489           density_y[y].direction[i].black+=
1490             cooccurrence[x][y].direction[i].black;
1491         if (image->alpha_trait == BlendPixelTrait)
1492           density_y[y].direction[i].alpha+=
1493             cooccurrence[x][y].direction[i].alpha;
1494       }
1495       mean.direction[i].red+=y*sum[y].direction[i].red;
1496       sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1497       mean.direction[i].green+=y*sum[y].direction[i].green;
1498       sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1499       mean.direction[i].blue+=y*sum[y].direction[i].blue;
1500       sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1501       if (image->colorspace == CMYKColorspace)
1502         {
1503           mean.direction[i].black+=y*sum[y].direction[i].black;
1504           sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1505         }
1506       if (image->alpha_trait == BlendPixelTrait)
1507         {
1508           mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1509           sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1510         }
1511     }
1512     /*
1513       Correlation: measure of linear-dependencies in the image.
1514     */
1515     channel_features[RedPixelChannel].correlation[i]=
1516       (correlation.direction[i].red-mean.direction[i].red*
1517       mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1518       (mean.direction[i].red*mean.direction[i].red))*sqrt(
1519       sum_squares.direction[i].red-(mean.direction[i].red*
1520       mean.direction[i].red)));
1521     channel_features[GreenPixelChannel].correlation[i]=
1522       (correlation.direction[i].green-mean.direction[i].green*
1523       mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1524       (mean.direction[i].green*mean.direction[i].green))*sqrt(
1525       sum_squares.direction[i].green-(mean.direction[i].green*
1526       mean.direction[i].green)));
1527     channel_features[BluePixelChannel].correlation[i]=
1528       (correlation.direction[i].blue-mean.direction[i].blue*
1529       mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1530       (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1531       sum_squares.direction[i].blue-(mean.direction[i].blue*
1532       mean.direction[i].blue)));
1533     if (image->colorspace == CMYKColorspace)
1534       channel_features[BlackPixelChannel].correlation[i]=
1535         (correlation.direction[i].black-mean.direction[i].black*
1536         mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1537         (mean.direction[i].black*mean.direction[i].black))*sqrt(
1538         sum_squares.direction[i].black-(mean.direction[i].black*
1539         mean.direction[i].black)));
1540     if (image->alpha_trait == BlendPixelTrait)
1541       channel_features[AlphaPixelChannel].correlation[i]=
1542         (correlation.direction[i].alpha-mean.direction[i].alpha*
1543         mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1544         (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1545         sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1546         mean.direction[i].alpha)));
1547   }
1548   /*
1549     Compute more texture features.
1550   */
1551 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1552   #pragma omp parallel for schedule(static,4) shared(status) \
1553     magick_threads(image,image,number_grays,1)
1554 #endif
1555   for (i=0; i < 4; i++)
1556   {
1557     register ssize_t
1558       x;
1559
1560     for (x=2; x < (ssize_t) (2*number_grays); x++)
1561     {
1562       /*
1563         Sum average.
1564       */
1565       channel_features[RedPixelChannel].sum_average[i]+=
1566         x*density_xy[x].direction[i].red;
1567       channel_features[GreenPixelChannel].sum_average[i]+=
1568         x*density_xy[x].direction[i].green;
1569       channel_features[BluePixelChannel].sum_average[i]+=
1570         x*density_xy[x].direction[i].blue;
1571       if (image->colorspace == CMYKColorspace)
1572         channel_features[BlackPixelChannel].sum_average[i]+=
1573           x*density_xy[x].direction[i].black;
1574       if (image->alpha_trait == BlendPixelTrait)
1575         channel_features[AlphaPixelChannel].sum_average[i]+=
1576           x*density_xy[x].direction[i].alpha;
1577       /*
1578         Sum entropy.
1579       */
1580       channel_features[RedPixelChannel].sum_entropy[i]-=
1581         density_xy[x].direction[i].red*
1582         MagickLog10(density_xy[x].direction[i].red);
1583       channel_features[GreenPixelChannel].sum_entropy[i]-=
1584         density_xy[x].direction[i].green*
1585         MagickLog10(density_xy[x].direction[i].green);
1586       channel_features[BluePixelChannel].sum_entropy[i]-=
1587         density_xy[x].direction[i].blue*
1588         MagickLog10(density_xy[x].direction[i].blue);
1589       if (image->colorspace == CMYKColorspace)
1590         channel_features[BlackPixelChannel].sum_entropy[i]-=
1591           density_xy[x].direction[i].black*
1592           MagickLog10(density_xy[x].direction[i].black);
1593       if (image->alpha_trait == BlendPixelTrait)
1594         channel_features[AlphaPixelChannel].sum_entropy[i]-=
1595           density_xy[x].direction[i].alpha*
1596           MagickLog10(density_xy[x].direction[i].alpha);
1597       /*
1598         Sum variance.
1599       */
1600       channel_features[RedPixelChannel].sum_variance[i]+=
1601         (x-channel_features[RedPixelChannel].sum_entropy[i])*
1602         (x-channel_features[RedPixelChannel].sum_entropy[i])*
1603         density_xy[x].direction[i].red;
1604       channel_features[GreenPixelChannel].sum_variance[i]+=
1605         (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1606         (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1607         density_xy[x].direction[i].green;
1608       channel_features[BluePixelChannel].sum_variance[i]+=
1609         (x-channel_features[BluePixelChannel].sum_entropy[i])*
1610         (x-channel_features[BluePixelChannel].sum_entropy[i])*
1611         density_xy[x].direction[i].blue;
1612       if (image->colorspace == CMYKColorspace)
1613         channel_features[BlackPixelChannel].sum_variance[i]+=
1614           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1615           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1616           density_xy[x].direction[i].black;
1617       if (image->alpha_trait == BlendPixelTrait)
1618         channel_features[AlphaPixelChannel].sum_variance[i]+=
1619           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1620           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1621           density_xy[x].direction[i].alpha;
1622     }
1623   }
1624   /*
1625     Compute more texture features.
1626   */
1627 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1628   #pragma omp parallel for schedule(static,4) shared(status) \
1629     magick_threads(image,image,number_grays,1)
1630 #endif
1631   for (i=0; i < 4; i++)
1632   {
1633     register ssize_t
1634       y;
1635
1636     for (y=0; y < (ssize_t) number_grays; y++)
1637     {
1638       register ssize_t
1639         x;
1640
1641       for (x=0; x < (ssize_t) number_grays; x++)
1642       {
1643         /*
1644           Sum of Squares: Variance
1645         */
1646         variance.direction[i].red+=(y-mean.direction[i].red+1)*
1647           (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1648         variance.direction[i].green+=(y-mean.direction[i].green+1)*
1649           (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1650         variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1651           (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1652         if (image->colorspace == CMYKColorspace)
1653           variance.direction[i].black+=(y-mean.direction[i].black+1)*
1654             (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1655         if (image->alpha_trait == BlendPixelTrait)
1656           variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1657             (y-mean.direction[i].alpha+1)*
1658             cooccurrence[x][y].direction[i].alpha;
1659         /*
1660           Sum average / Difference Variance.
1661         */
1662         density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1663           cooccurrence[x][y].direction[i].red;
1664         density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1665           cooccurrence[x][y].direction[i].green;
1666         density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1667           cooccurrence[x][y].direction[i].blue;
1668         if (image->colorspace == CMYKColorspace)
1669           density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1670             cooccurrence[x][y].direction[i].black;
1671         if (image->alpha_trait == BlendPixelTrait)
1672           density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1673             cooccurrence[x][y].direction[i].alpha;
1674         /*
1675           Information Measures of Correlation.
1676         */
1677         entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1678           MagickLog10(cooccurrence[x][y].direction[i].red);
1679         entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1680           MagickLog10(cooccurrence[x][y].direction[i].green);
1681         entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1682           MagickLog10(cooccurrence[x][y].direction[i].blue);
1683         if (image->colorspace == CMYKColorspace)
1684           entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1685             MagickLog10(cooccurrence[x][y].direction[i].black);
1686         if (image->alpha_trait == BlendPixelTrait)
1687           entropy_xy.direction[i].alpha-=
1688             cooccurrence[x][y].direction[i].alpha*MagickLog10(
1689             cooccurrence[x][y].direction[i].alpha);
1690         entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1691           MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
1692         entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1693           MagickLog10(density_x[x].direction[i].green*
1694           density_y[y].direction[i].green));
1695         entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1696           MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1697         if (image->colorspace == CMYKColorspace)
1698           entropy_xy1.direction[i].black-=(
1699             cooccurrence[x][y].direction[i].black*MagickLog10(
1700             density_x[x].direction[i].black*density_y[y].direction[i].black));
1701         if (image->alpha_trait == BlendPixelTrait)
1702           entropy_xy1.direction[i].alpha-=(
1703             cooccurrence[x][y].direction[i].alpha*MagickLog10(
1704             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1705         entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1706           density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
1707           density_y[y].direction[i].red));
1708         entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1709           density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
1710           density_y[y].direction[i].green));
1711         entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1712           density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
1713           density_y[y].direction[i].blue));
1714         if (image->colorspace == CMYKColorspace)
1715           entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1716             density_y[y].direction[i].black*MagickLog10(
1717             density_x[x].direction[i].black*density_y[y].direction[i].black));
1718         if (image->alpha_trait == BlendPixelTrait)
1719           entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1720             density_y[y].direction[i].alpha*MagickLog10(
1721             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1722       }
1723     }
1724     channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1725       variance.direction[i].red;
1726     channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1727       variance.direction[i].green;
1728     channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1729       variance.direction[i].blue;
1730     if (image->colorspace == CMYKColorspace)
1731       channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1732         variance.direction[i].black;
1733     if (image->alpha_trait == BlendPixelTrait)
1734       channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1735         variance.direction[i].alpha;
1736   }
1737   /*
1738     Compute more texture features.
1739   */
1740   (void) ResetMagickMemory(&variance,0,sizeof(variance));
1741   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
1742 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1743   #pragma omp parallel for schedule(static,4) shared(status) \
1744     magick_threads(image,image,number_grays,1)
1745 #endif
1746   for (i=0; i < 4; i++)
1747   {
1748     register ssize_t
1749       x;
1750
1751     for (x=0; x < (ssize_t) number_grays; x++)
1752     {
1753       /*
1754         Difference variance.
1755       */
1756       variance.direction[i].red+=density_xy[x].direction[i].red;
1757       variance.direction[i].green+=density_xy[x].direction[i].green;
1758       variance.direction[i].blue+=density_xy[x].direction[i].blue;
1759       if (image->colorspace == CMYKColorspace)
1760         variance.direction[i].black+=density_xy[x].direction[i].black;
1761       if (image->alpha_trait == BlendPixelTrait)
1762         variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1763       sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1764         density_xy[x].direction[i].red;
1765       sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1766         density_xy[x].direction[i].green;
1767       sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1768         density_xy[x].direction[i].blue;
1769       if (image->colorspace == CMYKColorspace)
1770         sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1771           density_xy[x].direction[i].black;
1772       if (image->alpha_trait == BlendPixelTrait)
1773         sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1774           density_xy[x].direction[i].alpha;
1775       /*
1776         Difference entropy.
1777       */
1778       channel_features[RedPixelChannel].difference_entropy[i]-=
1779         density_xy[x].direction[i].red*
1780         MagickLog10(density_xy[x].direction[i].red);
1781       channel_features[GreenPixelChannel].difference_entropy[i]-=
1782         density_xy[x].direction[i].green*
1783         MagickLog10(density_xy[x].direction[i].green);
1784       channel_features[BluePixelChannel].difference_entropy[i]-=
1785         density_xy[x].direction[i].blue*
1786         MagickLog10(density_xy[x].direction[i].blue);
1787       if (image->colorspace == CMYKColorspace)
1788         channel_features[BlackPixelChannel].difference_entropy[i]-=
1789           density_xy[x].direction[i].black*
1790           MagickLog10(density_xy[x].direction[i].black);
1791       if (image->alpha_trait == BlendPixelTrait)
1792         channel_features[AlphaPixelChannel].difference_entropy[i]-=
1793           density_xy[x].direction[i].alpha*
1794           MagickLog10(density_xy[x].direction[i].alpha);
1795       /*
1796         Information Measures of Correlation.
1797       */
1798       entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1799         MagickLog10(density_x[x].direction[i].red));
1800       entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1801         MagickLog10(density_x[x].direction[i].green));
1802       entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1803         MagickLog10(density_x[x].direction[i].blue));
1804       if (image->colorspace == CMYKColorspace)
1805         entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1806           MagickLog10(density_x[x].direction[i].black));
1807       if (image->alpha_trait == BlendPixelTrait)
1808         entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1809           MagickLog10(density_x[x].direction[i].alpha));
1810       entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1811         MagickLog10(density_y[x].direction[i].red));
1812       entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1813         MagickLog10(density_y[x].direction[i].green));
1814       entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1815         MagickLog10(density_y[x].direction[i].blue));
1816       if (image->colorspace == CMYKColorspace)
1817         entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1818           MagickLog10(density_y[x].direction[i].black));
1819       if (image->alpha_trait == BlendPixelTrait)
1820         entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1821           MagickLog10(density_y[x].direction[i].alpha));
1822     }
1823     /*
1824       Difference variance.
1825     */
1826     channel_features[RedPixelChannel].difference_variance[i]=
1827       (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1828       (variance.direction[i].red*variance.direction[i].red))/
1829       ((double) number_grays*number_grays*number_grays*number_grays);
1830     channel_features[GreenPixelChannel].difference_variance[i]=
1831       (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1832       (variance.direction[i].green*variance.direction[i].green))/
1833       ((double) number_grays*number_grays*number_grays*number_grays);
1834     channel_features[BluePixelChannel].difference_variance[i]=
1835       (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1836       (variance.direction[i].blue*variance.direction[i].blue))/
1837       ((double) number_grays*number_grays*number_grays*number_grays);
1838     if (image->colorspace == CMYKColorspace)
1839       channel_features[BlackPixelChannel].difference_variance[i]=
1840         (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1841         (variance.direction[i].black*variance.direction[i].black))/
1842         ((double) number_grays*number_grays*number_grays*number_grays);
1843     if (image->alpha_trait == BlendPixelTrait)
1844       channel_features[AlphaPixelChannel].difference_variance[i]=
1845         (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1846         (variance.direction[i].alpha*variance.direction[i].alpha))/
1847         ((double) number_grays*number_grays*number_grays*number_grays);
1848     /*
1849       Information Measures of Correlation.
1850     */
1851     channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1852       (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1853       (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1854        entropy_x.direction[i].red : entropy_y.direction[i].red);
1855     channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1856       (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1857       (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1858        entropy_x.direction[i].green : entropy_y.direction[i].green);
1859     channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1860       (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1861       (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1862        entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1863     if (image->colorspace == CMYKColorspace)
1864       channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1865         (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1866         (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1867          entropy_x.direction[i].black : entropy_y.direction[i].black);
1868     if (image->alpha_trait == BlendPixelTrait)
1869       channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1870         (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1871         (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1872          entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1873     channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1874       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
1875       entropy_xy.direction[i].red)))));
1876     channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1877       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
1878       entropy_xy.direction[i].green)))));
1879     channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1880       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
1881       entropy_xy.direction[i].blue)))));
1882     if (image->colorspace == CMYKColorspace)
1883       channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1884         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black-
1885         entropy_xy.direction[i].black)))));
1886     if (image->alpha_trait == BlendPixelTrait)
1887       channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1888         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha-
1889         entropy_xy.direction[i].alpha)))));
1890   }
1891   /*
1892     Compute more texture features.
1893   */
1894 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1895   #pragma omp parallel for schedule(static,4) shared(status) \
1896     magick_threads(image,image,number_grays,1)
1897 #endif
1898   for (i=0; i < 4; i++)
1899   {
1900     ssize_t
1901       z;
1902
1903     for (z=0; z < (ssize_t) number_grays; z++)
1904     {
1905       register ssize_t
1906         y;
1907
1908       ChannelStatistics
1909         pixel;
1910
1911       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
1912       for (y=0; y < (ssize_t) number_grays; y++)
1913       {
1914         register ssize_t
1915           x;
1916
1917         for (x=0; x < (ssize_t) number_grays; x++)
1918         {
1919           /*
1920             Contrast:  amount of local variations present in an image.
1921           */
1922           if (((y-x) == z) || ((x-y) == z))
1923             {
1924               pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1925               pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1926               pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1927               if (image->colorspace == CMYKColorspace)
1928                 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1929               if (image->alpha_trait == BlendPixelTrait)
1930                 pixel.direction[i].alpha+=
1931                   cooccurrence[x][y].direction[i].alpha;
1932             }
1933           /*
1934             Maximum Correlation Coefficient.
1935           */
1936           Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1937             cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1938             density_y[x].direction[i].red;
1939           Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1940             cooccurrence[y][x].direction[i].green/
1941             density_x[z].direction[i].green/density_y[x].direction[i].red;
1942           Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1943             cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
1944             density_y[x].direction[i].blue;
1945           if (image->colorspace == CMYKColorspace)
1946             Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1947               cooccurrence[y][x].direction[i].black/
1948               density_x[z].direction[i].black/density_y[x].direction[i].black;
1949           if (image->alpha_trait == BlendPixelTrait)
1950             Q[z][y].direction[i].alpha+=
1951               cooccurrence[z][x].direction[i].alpha*
1952               cooccurrence[y][x].direction[i].alpha/
1953               density_x[z].direction[i].alpha/
1954               density_y[x].direction[i].alpha;
1955         }
1956       }
1957       channel_features[RedPixelChannel].contrast[i]+=z*z*
1958         pixel.direction[i].red;
1959       channel_features[GreenPixelChannel].contrast[i]+=z*z*
1960         pixel.direction[i].green;
1961       channel_features[BluePixelChannel].contrast[i]+=z*z*
1962         pixel.direction[i].blue;
1963       if (image->colorspace == CMYKColorspace)
1964         channel_features[BlackPixelChannel].contrast[i]+=z*z*
1965           pixel.direction[i].black;
1966       if (image->alpha_trait == BlendPixelTrait)
1967         channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1968           pixel.direction[i].alpha;
1969     }
1970     /*
1971       Maximum Correlation Coefficient.
1972       Future: return second largest eigenvalue of Q.
1973     */
1974     channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1975       sqrt((double) -1.0);
1976     channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1977       sqrt((double) -1.0);
1978     channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1979       sqrt((double) -1.0);
1980     if (image->colorspace == CMYKColorspace)
1981       channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1982         sqrt((double) -1.0);
1983     if (image->alpha_trait == BlendPixelTrait)
1984       channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1985         sqrt((double) -1.0);
1986   }
1987   /*
1988     Relinquish resources.
1989   */
1990   sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1991   for (i=0; i < (ssize_t) number_grays; i++)
1992     Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1993   Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1994   density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1995   density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1996   density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1997   for (i=0; i < (ssize_t) number_grays; i++)
1998     cooccurrence[i]=(ChannelStatistics *)
1999       RelinquishMagickMemory(cooccurrence[i]);
2000   cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
2001   return(channel_features);
2002 }