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,"viewbox 0 0 %.20g %.20g\n",
 722     (double) image->columns,(double) image->rows);
 723   (void) write(file,message,strlen(message));
 724   line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
 725   if (threshold != 0)
 726     line_count=threshold;
 727   for (y=0; y < (ssize_t) accumulator_height; y++)
 728   {
 729     register ssize_t
 730       x;
 731
 732     for (x=0; x < (ssize_t) accumulator_width; x++)
 733     {
 734       double
 735         count;
 736
 737       (void) GetMatrixElement(accumulator,x,y,&count);
 738       if (count >= (double) line_count)
 739         {
 740           double
 741             maxima,
 742             x1,
 743             x2,
 744             y1,
 745             y2;
 746
 747           ssize_t
 748             v;
 749
 750           /*
 751             Is point a local maxima?
 752           */
 753           maxima=count;
 754           for (v=((ssize_t) -(height/2)); v < ((ssize_t) (height/2)); v++)
 755           {
 756             ssize_t
 757               u;
 758
 759             for (u=((ssize_t) -(width/2)); u < ((ssize_t) (width/2)); u++)
 760             {
 761               if ((u != 0) || (v !=0))
 762                 {
 763                   (void) GetMatrixElement(accumulator,x+u,y+v,&count);
 764                   if (count > maxima)
 765                     {
 766                       maxima=count;
 767                       break;
 768                     }
 769                 }
 770             }
 771             if (u < (ssize_t) (width/2))
 772               break;
 773           }
 774           (void) GetMatrixElement(accumulator,x,y,&count);
 775           if (maxima > count)
 776             continue;
 777           if ((x >= 45) && (x <= 135))
 778             {
 779               /*
 780                 y = (r-x cos(t))/sin(t)
 781               */
 782               x1=0.0;
 783               y1=(double) ((y-(accumulator_height/2.0))-((x1-(image->columns/
 784                 2.0))*cos(DegreesToRadians((double) x))))/
 785                 sin(DegreesToRadians((double) x))+(image->rows/2.0);
 786               x2=(double) image->columns;
 787               y2=(double) ((y-(accumulator_height/2.0))-((x2-(image->columns/
 788                 2.0))*cos(DegreesToRadians((double) x))))/
 789                 sin(DegreesToRadians((double) x))+(image->rows/2.0);
 790             }
 791           else
 792             {
 793               /*
 794                 x = (r-y cos(t))/sin(t)
 795               */
 796               y1=0.0;
 797               x1=(double) ((y-(accumulator_height/2.0))-((y1-(image->rows/2.0))*
 798                 sin(DegreesToRadians((double) x))))/cos(DegreesToRadians(
 799                 (double) x))+(image->columns/2.0);
 800               y2=(double) image->rows;
 801               x2=(double) ((y-(accumulator_height/2.0))-((y2-(image->rows/2.0))*
 802                 sin(DegreesToRadians((double) x))))/cos(DegreesToRadians(
 803                 (double) x))+(image->columns/2.0);
 804             }
 805           (void) FormatLocaleString(message,MaxTextExtent,"line %g,%g %g,%g\n",
 806             x1,y1,x2,y2);
 807           (void) write(file,message,strlen(message));
 808         }
 809     }
 810   }
 811   (void) close(file);
 812   /*
 813     Render lines to image canvas.
 814   */
 815   image_info=AcquireImageInfo();
 816   image_info->background_color=image->background_color;
 817   (void) FormatLocaleString(image_info->filename,MaxTextExtent,"mvg:%s",path);
 818   artifact=GetImageArtifact(image,"background");
 819   if (artifact != (const char *) NULL)
 820     (void) SetImageOption(image_info,"background",artifact);
 821   artifact=GetImageArtifact(image,"fill");
 822   if (artifact != (const char *) NULL)
 823     (void) SetImageOption(image_info,"fill",artifact);
 824   artifact=GetImageArtifact(image,"stroke");
 825   if (artifact != (const char *) NULL)
 826     (void) SetImageOption(image_info,"stroke",artifact);
 827   artifact=GetImageArtifact(image,"strokewidth");
 828   if (artifact != (const char *) NULL)
 829     (void) SetImageOption(image_info,"strokewidth",artifact);
 830   lines_image=ReadImage(image_info,exception);
 831   artifact=GetImageArtifact(image,"hough-lines:accumulator");
 832   if ((lines_image != (Image *) NULL) &&
 833       (IsStringTrue(artifact) != MagickFalse))
 834     {
 835       Image
 836         *accumulator_image;
 837
 838       accumulator_image=MatrixToImage(accumulator,exception);
 839       if (accumulator_image != (Image *) NULL)
 840         AppendImageToList(&lines_image,accumulator_image);
 841     }
 842   /*
 843     Free resources.
 844   */
 845   accumulator=DestroyMatrixInfo(accumulator);
 846   image_info=DestroyImageInfo(image_info);
 847   (void) RelinquishUniqueFileResource(path);
 848   return(GetFirstImageInList(lines_image));
 849 }
 850 \f
 851 /*
 852 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 853 %                                                                             %
 854 %                                                                             %
 855 %                                                                             %
 856 %   G e t I m a g e F e a t u r e s                                           %
 857 %                                                                             %
 858 %                                                                             %
 859 %                                                                             %
 860 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 861 %
 862 %  GetImageFeatures() returns features for each channel in the image in
 863 %  each of four directions (horizontal, vertical, left and right diagonals)
 864 %  for the specified distance.  The features include the angular second
 865 %  moment, contrast, correlation, sum of squares: variance, inverse difference
 866 %  moment, sum average, sum varience, sum entropy, entropy, difference variance,%  difference entropy, information measures of correlation 1, information
 867 %  measures of correlation 2, and maximum correlation coefficient.  You can
 868 %  access the red channel contrast, for example, like this:
 869 %
 870 %      channel_features=GetImageFeatures(image,1,exception);
 871 %      contrast=channel_features[RedPixelChannel].contrast[0];
 872 %
 873 %  Use MagickRelinquishMemory() to free the features buffer.
 874 %
 875 %  The format of the GetImageFeatures method is:
 876 %
 877 %      ChannelFeatures *GetImageFeatures(const Image *image,
 878 %        const size_t distance,ExceptionInfo *exception)
 879 %
 880 %  A description of each parameter follows:
 881 %
 882 %    o image: the image.
 883 %
 884 %    o distance: the distance.
 885 %
 886 %    o exception: return any errors or warnings in this structure.
 887 %
 888 */
 889
 890 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
 891 {
 892   if (x < 0)
 893     return(-x);
 894   return(x);
 895 }
 896
 897 static inline double MagickLog10(const double x)
 898 {
 899 #define Log10Epsilon  (1.0e-11)
 900
 901  if (fabs(x) < Log10Epsilon)
 902    return(log10(Log10Epsilon));
 903  return(log10(fabs(x)));
 904 }
 905
 906 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
 907   const size_t distance,ExceptionInfo *exception)
 908 {
 909   typedef struct _ChannelStatistics
 910   {
 911     PixelInfo
 912       direction[4];  /* horizontal, vertical, left and right diagonals */
 913   } ChannelStatistics;
 914
 915   CacheView
 916     *image_view;
 917
 918   ChannelFeatures
 919     *channel_features;
 920
 921   ChannelStatistics
 922     **cooccurrence,
 923     correlation,
 924     *density_x,
 925     *density_xy,
 926     *density_y,
 927     entropy_x,
 928     entropy_xy,
 929     entropy_xy1,
 930     entropy_xy2,
 931     entropy_y,
 932     mean,
 933     **Q,
 934     *sum,
 935     sum_squares,
 936     variance;
 937
 938   PixelPacket
 939     gray,
 940     *grays;
 941
 942   MagickBooleanType
 943     status;
 944
 945   register ssize_t
 946     i;
 947
 948   size_t
 949     length;
 950
 951   ssize_t
 952     y;
 953
 954   unsigned int
 955     number_grays;
 956
 957   assert(image != (Image *) NULL);
 958   assert(image->signature == MagickSignature);
 959   if (image->debug != MagickFalse)
 960     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 961   if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
 962     return((ChannelFeatures *) NULL);
 963   length=CompositeChannels+1UL;
 964   channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
 965     sizeof(*channel_features));
 966   if (channel_features == (ChannelFeatures *) NULL)
 967     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 968   (void) ResetMagickMemory(channel_features,0,length*
 969     sizeof(*channel_features));
 970   /*
 971     Form grays.
 972   */
 973   grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
 974   if (grays == (PixelPacket *) NULL)
 975     {
 976       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 977         channel_features);
 978       (void) ThrowMagickException(exception,GetMagickModule(),
 979         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 980       return(channel_features);
 981     }
 982   for (i=0; i <= (ssize_t) MaxMap; i++)
 983   {
 984     grays[i].red=(~0U);
 985     grays[i].green=(~0U);
 986     grays[i].blue=(~0U);
 987     grays[i].alpha=(~0U);
 988     grays[i].black=(~0U);
 989   }
 990   status=MagickTrue;
 991   image_view=AcquireVirtualCacheView(image,exception);
 992 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 993   #pragma omp parallel for schedule(static,4) shared(status) \
 994     magick_threads(image,image,image->rows,1)
 995 #endif
 996   for (y=0; y < (ssize_t) image->rows; y++)
 997   {
 998     register const Quantum
 999       *restrict p;
1000
1001     register ssize_t
1002       x;
1003
1004     if (status == MagickFalse)
1005       continue;
1006     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1007     if (p == (const Quantum *) NULL)
1008       {
1009         status=MagickFalse;
1010         continue;
1011       }
1012     for (x=0; x < (ssize_t) image->columns; x++)
1013     {
1014       grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
1015         ScaleQuantumToMap(GetPixelRed(image,p));
1016       grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
1017         ScaleQuantumToMap(GetPixelGreen(image,p));
1018       grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
1019         ScaleQuantumToMap(GetPixelBlue(image,p));
1020       if (image->colorspace == CMYKColorspace)
1021         grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
1022           ScaleQuantumToMap(GetPixelBlack(image,p));
1023       if (image->alpha_trait == BlendPixelTrait)
1024         grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
1025           ScaleQuantumToMap(GetPixelAlpha(image,p));
1026       p+=GetPixelChannels(image);
1027     }
1028   }
1029   image_view=DestroyCacheView(image_view);
1030   if (status == MagickFalse)
1031     {
1032       grays=(PixelPacket *) RelinquishMagickMemory(grays);
1033       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1034         channel_features);
1035       return(channel_features);
1036     }
1037   (void) ResetMagickMemory(&gray,0,sizeof(gray));
1038   for (i=0; i <= (ssize_t) MaxMap; i++)
1039   {
1040     if (grays[i].red != ~0U)
1041       grays[gray.red++].red=grays[i].red;
1042     if (grays[i].green != ~0U)
1043       grays[gray.green++].green=grays[i].green;
1044     if (grays[i].blue != ~0U)
1045       grays[gray.blue++].blue=grays[i].blue;
1046     if (image->colorspace == CMYKColorspace)
1047       if (grays[i].black != ~0U)
1048         grays[gray.black++].black=grays[i].black;
1049     if (image->alpha_trait == BlendPixelTrait)
1050       if (grays[i].alpha != ~0U)
1051         grays[gray.alpha++].alpha=grays[i].alpha;
1052   }
1053   /*
1054     Allocate spatial dependence matrix.
1055   */
1056   number_grays=gray.red;
1057   if (gray.green > number_grays)
1058     number_grays=gray.green;
1059   if (gray.blue > number_grays)
1060     number_grays=gray.blue;
1061   if (image->colorspace == CMYKColorspace)
1062     if (gray.black > number_grays)
1063       number_grays=gray.black;
1064   if (image->alpha_trait == BlendPixelTrait)
1065     if (gray.alpha > number_grays)
1066       number_grays=gray.alpha;
1067   cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
1068     sizeof(*cooccurrence));
1069   density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
1070     sizeof(*density_x));
1071   density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
1072     sizeof(*density_xy));
1073   density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
1074     sizeof(*density_y));
1075   Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
1076   sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
1077   if ((cooccurrence == (ChannelStatistics **) NULL) ||
1078       (density_x == (ChannelStatistics *) NULL) ||
1079       (density_xy == (ChannelStatistics *) NULL) ||
1080       (density_y == (ChannelStatistics *) NULL) ||
1081       (Q == (ChannelStatistics **) NULL) ||
1082       (sum == (ChannelStatistics *) NULL))
1083     {
1084       if (Q != (ChannelStatistics **) NULL)
1085         {
1086           for (i=0; i < (ssize_t) number_grays; i++)
1087             Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1088           Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1089         }
1090       if (sum != (ChannelStatistics *) NULL)
1091         sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1092       if (density_y != (ChannelStatistics *) NULL)
1093         density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1094       if (density_xy != (ChannelStatistics *) NULL)
1095         density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1096       if (density_x != (ChannelStatistics *) NULL)
1097         density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1098       if (cooccurrence != (ChannelStatistics **) NULL)
1099         {
1100           for (i=0; i < (ssize_t) number_grays; i++)
1101             cooccurrence[i]=(ChannelStatistics *)
1102               RelinquishMagickMemory(cooccurrence[i]);
1103           cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
1104             cooccurrence);
1105         }
1106       grays=(PixelPacket *) RelinquishMagickMemory(grays);
1107       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1108         channel_features);
1109       (void) ThrowMagickException(exception,GetMagickModule(),
1110         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1111       return(channel_features);
1112     }
1113   (void) ResetMagickMemory(&correlation,0,sizeof(correlation));
1114   (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x));
1115   (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
1116   (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y));
1117   (void) ResetMagickMemory(&mean,0,sizeof(mean));
1118   (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum));
1119   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
1120   (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy));
1121   (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x));
1122   (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy));
1123   (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1));
1124   (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2));
1125   (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y));
1126   (void) ResetMagickMemory(&variance,0,sizeof(variance));
1127   for (i=0; i < (ssize_t) number_grays; i++)
1128   {
1129     cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
1130       sizeof(**cooccurrence));
1131     Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
1132     if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
1133         (Q[i] == (ChannelStatistics *) NULL))
1134       break;
1135     (void) ResetMagickMemory(cooccurrence[i],0,number_grays*
1136       sizeof(**cooccurrence));
1137     (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q));
1138   }
1139   if (i < (ssize_t) number_grays)
1140     {
1141       for (i--; i >= 0; i--)
1142       {
1143         if (Q[i] != (ChannelStatistics *) NULL)
1144           Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1145         if (cooccurrence[i] != (ChannelStatistics *) NULL)
1146           cooccurrence[i]=(ChannelStatistics *)
1147             RelinquishMagickMemory(cooccurrence[i]);
1148       }
1149       Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1150       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1151       sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1152       density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1153       density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1154       density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1155       grays=(PixelPacket *) RelinquishMagickMemory(grays);
1156       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1157         channel_features);
1158       (void) ThrowMagickException(exception,GetMagickModule(),
1159         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1160       return(channel_features);
1161     }
1162   /*
1163     Initialize spatial dependence matrix.
1164   */
1165   status=MagickTrue;
1166   image_view=AcquireVirtualCacheView(image,exception);
1167   for (y=0; y < (ssize_t) image->rows; y++)
1168   {
1169     register const Quantum
1170       *restrict p;
1171
1172     register ssize_t
1173       x;
1174
1175     ssize_t
1176       i,
1177       offset,
1178       u,
1179       v;
1180
1181     if (status == MagickFalse)
1182       continue;
1183     p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+
1184       2*distance,distance+2,exception);
1185     if (p == (const Quantum *) NULL)
1186       {
1187         status=MagickFalse;
1188         continue;
1189       }
1190     p+=distance*GetPixelChannels(image);;
1191     for (x=0; x < (ssize_t) image->columns; x++)
1192     {
1193       for (i=0; i < 4; i++)
1194       {
1195         switch (i)
1196         {
1197           case 0:
1198           default:
1199           {
1200             /*
1201               Horizontal adjacency.
1202             */
1203             offset=(ssize_t) distance;
1204             break;
1205           }
1206           case 1:
1207           {
1208             /*
1209               Vertical adjacency.
1210             */
1211             offset=(ssize_t) (image->columns+2*distance);
1212             break;
1213           }
1214           case 2:
1215           {
1216             /*
1217               Right diagonal adjacency.
1218             */
1219             offset=(ssize_t) ((image->columns+2*distance)-distance);
1220             break;
1221           }
1222           case 3:
1223           {
1224             /*
1225               Left diagonal adjacency.
1226             */
1227             offset=(ssize_t) ((image->columns+2*distance)+distance);
1228             break;
1229           }
1230         }
1231         u=0;
1232         v=0;
1233         while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
1234           u++;
1235         while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
1236           v++;
1237         cooccurrence[u][v].direction[i].red++;
1238         cooccurrence[v][u].direction[i].red++;
1239         u=0;
1240         v=0;
1241         while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
1242           u++;
1243         while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
1244           v++;
1245         cooccurrence[u][v].direction[i].green++;
1246         cooccurrence[v][u].direction[i].green++;
1247         u=0;
1248         v=0;
1249         while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
1250           u++;
1251         while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
1252           v++;
1253         cooccurrence[u][v].direction[i].blue++;
1254         cooccurrence[v][u].direction[i].blue++;
1255         if (image->colorspace == CMYKColorspace)
1256           {
1257             u=0;
1258             v=0;
1259             while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
1260               u++;
1261             while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
1262               v++;
1263             cooccurrence[u][v].direction[i].black++;
1264             cooccurrence[v][u].direction[i].black++;
1265           }
1266         if (image->alpha_trait == BlendPixelTrait)
1267           {
1268             u=0;
1269             v=0;
1270             while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
1271               u++;
1272             while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
1273               v++;
1274             cooccurrence[u][v].direction[i].alpha++;
1275             cooccurrence[v][u].direction[i].alpha++;
1276           }
1277       }
1278       p+=GetPixelChannels(image);
1279     }
1280   }
1281   grays=(PixelPacket *) RelinquishMagickMemory(grays);
1282   image_view=DestroyCacheView(image_view);
1283   if (status == MagickFalse)
1284     {
1285       for (i=0; i < (ssize_t) number_grays; i++)
1286         cooccurrence[i]=(ChannelStatistics *)
1287           RelinquishMagickMemory(cooccurrence[i]);
1288       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1289       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1290         channel_features);
1291       (void) ThrowMagickException(exception,GetMagickModule(),
1292         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1293       return(channel_features);
1294     }
1295   /*
1296     Normalize spatial dependence matrix.
1297   */
1298   for (i=0; i < 4; i++)
1299   {
1300     double
1301       normalize;
1302
1303     register ssize_t
1304       y;
1305
1306     switch (i)
1307     {
1308       case 0:
1309       default:
1310       {
1311         /*
1312           Horizontal adjacency.
1313         */
1314         normalize=2.0*image->rows*(image->columns-distance);
1315         break;
1316       }
1317       case 1:
1318       {
1319         /*
1320           Vertical adjacency.
1321         */
1322         normalize=2.0*(image->rows-distance)*image->columns;
1323         break;
1324       }
1325       case 2:
1326       {
1327         /*
1328           Right diagonal adjacency.
1329         */
1330         normalize=2.0*(image->rows-distance)*(image->columns-distance);
1331         break;
1332       }
1333       case 3:
1334       {
1335         /*
1336           Left diagonal adjacency.
1337         */
1338         normalize=2.0*(image->rows-distance)*(image->columns-distance);
1339         break;
1340       }
1341     }
1342     normalize=PerceptibleReciprocal(normalize);
1343     for (y=0; y < (ssize_t) number_grays; y++)
1344     {
1345       register ssize_t
1346         x;
1347
1348       for (x=0; x < (ssize_t) number_grays; x++)
1349       {
1350         cooccurrence[x][y].direction[i].red*=normalize;
1351         cooccurrence[x][y].direction[i].green*=normalize;
1352         cooccurrence[x][y].direction[i].blue*=normalize;
1353         if (image->colorspace == CMYKColorspace)
1354           cooccurrence[x][y].direction[i].black*=normalize;
1355         if (image->alpha_trait == BlendPixelTrait)
1356           cooccurrence[x][y].direction[i].alpha*=normalize;
1357       }
1358     }
1359   }
1360   /*
1361     Compute texture features.
1362   */
1363 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1364   #pragma omp parallel for schedule(static,4) shared(status) \
1365     magick_threads(image,image,number_grays,1)
1366 #endif
1367   for (i=0; i < 4; i++)
1368   {
1369     register ssize_t
1370       y;
1371
1372     for (y=0; y < (ssize_t) number_grays; y++)
1373     {
1374       register ssize_t
1375         x;
1376
1377       for (x=0; x < (ssize_t) number_grays; x++)
1378       {
1379         /*
1380           Angular second moment:  measure of homogeneity of the image.
1381         */
1382         channel_features[RedPixelChannel].angular_second_moment[i]+=
1383           cooccurrence[x][y].direction[i].red*
1384           cooccurrence[x][y].direction[i].red;
1385         channel_features[GreenPixelChannel].angular_second_moment[i]+=
1386           cooccurrence[x][y].direction[i].green*
1387           cooccurrence[x][y].direction[i].green;
1388         channel_features[BluePixelChannel].angular_second_moment[i]+=
1389           cooccurrence[x][y].direction[i].blue*
1390           cooccurrence[x][y].direction[i].blue;
1391         if (image->colorspace == CMYKColorspace)
1392           channel_features[BlackPixelChannel].angular_second_moment[i]+=
1393             cooccurrence[x][y].direction[i].black*
1394             cooccurrence[x][y].direction[i].black;
1395         if (image->alpha_trait == BlendPixelTrait)
1396           channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1397             cooccurrence[x][y].direction[i].alpha*
1398             cooccurrence[x][y].direction[i].alpha;
1399         /*
1400           Correlation: measure of linear-dependencies in the image.
1401         */
1402         sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1403         sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1404         sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1405         if (image->colorspace == CMYKColorspace)
1406           sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1407         if (image->alpha_trait == BlendPixelTrait)
1408           sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1409         correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1410         correlation.direction[i].green+=x*y*
1411           cooccurrence[x][y].direction[i].green;
1412         correlation.direction[i].blue+=x*y*
1413           cooccurrence[x][y].direction[i].blue;
1414         if (image->colorspace == CMYKColorspace)
1415           correlation.direction[i].black+=x*y*
1416             cooccurrence[x][y].direction[i].black;
1417         if (image->alpha_trait == BlendPixelTrait)
1418           correlation.direction[i].alpha+=x*y*
1419             cooccurrence[x][y].direction[i].alpha;
1420         /*
1421           Inverse Difference Moment.
1422         */
1423         channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1424           cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1425         channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1426           cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1427         channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1428           cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1429         if (image->colorspace == CMYKColorspace)
1430           channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1431             cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1432         if (image->alpha_trait == BlendPixelTrait)
1433           channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1434             cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1435         /*
1436           Sum average.
1437         */
1438         density_xy[y+x+2].direction[i].red+=
1439           cooccurrence[x][y].direction[i].red;
1440         density_xy[y+x+2].direction[i].green+=
1441           cooccurrence[x][y].direction[i].green;
1442         density_xy[y+x+2].direction[i].blue+=
1443           cooccurrence[x][y].direction[i].blue;
1444         if (image->colorspace == CMYKColorspace)
1445           density_xy[y+x+2].direction[i].black+=
1446             cooccurrence[x][y].direction[i].black;
1447         if (image->alpha_trait == BlendPixelTrait)
1448           density_xy[y+x+2].direction[i].alpha+=
1449             cooccurrence[x][y].direction[i].alpha;
1450         /*
1451           Entropy.
1452         */
1453         channel_features[RedPixelChannel].entropy[i]-=
1454           cooccurrence[x][y].direction[i].red*
1455           MagickLog10(cooccurrence[x][y].direction[i].red);
1456         channel_features[GreenPixelChannel].entropy[i]-=
1457           cooccurrence[x][y].direction[i].green*
1458           MagickLog10(cooccurrence[x][y].direction[i].green);
1459         channel_features[BluePixelChannel].entropy[i]-=
1460           cooccurrence[x][y].direction[i].blue*
1461           MagickLog10(cooccurrence[x][y].direction[i].blue);
1462         if (image->colorspace == CMYKColorspace)
1463           channel_features[BlackPixelChannel].entropy[i]-=
1464             cooccurrence[x][y].direction[i].black*
1465             MagickLog10(cooccurrence[x][y].direction[i].black);
1466         if (image->alpha_trait == BlendPixelTrait)
1467           channel_features[AlphaPixelChannel].entropy[i]-=
1468             cooccurrence[x][y].direction[i].alpha*
1469             MagickLog10(cooccurrence[x][y].direction[i].alpha);
1470         /*
1471           Information Measures of Correlation.
1472         */
1473         density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1474         density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1475         density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1476         if (image->alpha_trait == BlendPixelTrait)
1477           density_x[x].direction[i].alpha+=
1478             cooccurrence[x][y].direction[i].alpha;
1479         if (image->colorspace == CMYKColorspace)
1480           density_x[x].direction[i].black+=
1481             cooccurrence[x][y].direction[i].black;
1482         density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1483         density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1484         density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1485         if (image->colorspace == CMYKColorspace)
1486           density_y[y].direction[i].black+=
1487             cooccurrence[x][y].direction[i].black;
1488         if (image->alpha_trait == BlendPixelTrait)
1489           density_y[y].direction[i].alpha+=
1490             cooccurrence[x][y].direction[i].alpha;
1491       }
1492       mean.direction[i].red+=y*sum[y].direction[i].red;
1493       sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1494       mean.direction[i].green+=y*sum[y].direction[i].green;
1495       sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1496       mean.direction[i].blue+=y*sum[y].direction[i].blue;
1497       sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1498       if (image->colorspace == CMYKColorspace)
1499         {
1500           mean.direction[i].black+=y*sum[y].direction[i].black;
1501           sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1502         }
1503       if (image->alpha_trait == BlendPixelTrait)
1504         {
1505           mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1506           sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1507         }
1508     }
1509     /*
1510       Correlation: measure of linear-dependencies in the image.
1511     */
1512     channel_features[RedPixelChannel].correlation[i]=
1513       (correlation.direction[i].red-mean.direction[i].red*
1514       mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1515       (mean.direction[i].red*mean.direction[i].red))*sqrt(
1516       sum_squares.direction[i].red-(mean.direction[i].red*
1517       mean.direction[i].red)));
1518     channel_features[GreenPixelChannel].correlation[i]=
1519       (correlation.direction[i].green-mean.direction[i].green*
1520       mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1521       (mean.direction[i].green*mean.direction[i].green))*sqrt(
1522       sum_squares.direction[i].green-(mean.direction[i].green*
1523       mean.direction[i].green)));
1524     channel_features[BluePixelChannel].correlation[i]=
1525       (correlation.direction[i].blue-mean.direction[i].blue*
1526       mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1527       (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1528       sum_squares.direction[i].blue-(mean.direction[i].blue*
1529       mean.direction[i].blue)));
1530     if (image->colorspace == CMYKColorspace)
1531       channel_features[BlackPixelChannel].correlation[i]=
1532         (correlation.direction[i].black-mean.direction[i].black*
1533         mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1534         (mean.direction[i].black*mean.direction[i].black))*sqrt(
1535         sum_squares.direction[i].black-(mean.direction[i].black*
1536         mean.direction[i].black)));
1537     if (image->alpha_trait == BlendPixelTrait)
1538       channel_features[AlphaPixelChannel].correlation[i]=
1539         (correlation.direction[i].alpha-mean.direction[i].alpha*
1540         mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1541         (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1542         sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1543         mean.direction[i].alpha)));
1544   }
1545   /*
1546     Compute more texture features.
1547   */
1548 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1549   #pragma omp parallel for schedule(static,4) shared(status) \
1550     magick_threads(image,image,number_grays,1)
1551 #endif
1552   for (i=0; i < 4; i++)
1553   {
1554     register ssize_t
1555       x;
1556
1557     for (x=2; x < (ssize_t) (2*number_grays); x++)
1558     {
1559       /*
1560         Sum average.
1561       */
1562       channel_features[RedPixelChannel].sum_average[i]+=
1563         x*density_xy[x].direction[i].red;
1564       channel_features[GreenPixelChannel].sum_average[i]+=
1565         x*density_xy[x].direction[i].green;
1566       channel_features[BluePixelChannel].sum_average[i]+=
1567         x*density_xy[x].direction[i].blue;
1568       if (image->colorspace == CMYKColorspace)
1569         channel_features[BlackPixelChannel].sum_average[i]+=
1570           x*density_xy[x].direction[i].black;
1571       if (image->alpha_trait == BlendPixelTrait)
1572         channel_features[AlphaPixelChannel].sum_average[i]+=
1573           x*density_xy[x].direction[i].alpha;
1574       /*
1575         Sum entropy.
1576       */
1577       channel_features[RedPixelChannel].sum_entropy[i]-=
1578         density_xy[x].direction[i].red*
1579         MagickLog10(density_xy[x].direction[i].red);
1580       channel_features[GreenPixelChannel].sum_entropy[i]-=
1581         density_xy[x].direction[i].green*
1582         MagickLog10(density_xy[x].direction[i].green);
1583       channel_features[BluePixelChannel].sum_entropy[i]-=
1584         density_xy[x].direction[i].blue*
1585         MagickLog10(density_xy[x].direction[i].blue);
1586       if (image->colorspace == CMYKColorspace)
1587         channel_features[BlackPixelChannel].sum_entropy[i]-=
1588           density_xy[x].direction[i].black*
1589           MagickLog10(density_xy[x].direction[i].black);
1590       if (image->alpha_trait == BlendPixelTrait)
1591         channel_features[AlphaPixelChannel].sum_entropy[i]-=
1592           density_xy[x].direction[i].alpha*
1593           MagickLog10(density_xy[x].direction[i].alpha);
1594       /*
1595         Sum variance.
1596       */
1597       channel_features[RedPixelChannel].sum_variance[i]+=
1598         (x-channel_features[RedPixelChannel].sum_entropy[i])*
1599         (x-channel_features[RedPixelChannel].sum_entropy[i])*
1600         density_xy[x].direction[i].red;
1601       channel_features[GreenPixelChannel].sum_variance[i]+=
1602         (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1603         (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1604         density_xy[x].direction[i].green;
1605       channel_features[BluePixelChannel].sum_variance[i]+=
1606         (x-channel_features[BluePixelChannel].sum_entropy[i])*
1607         (x-channel_features[BluePixelChannel].sum_entropy[i])*
1608         density_xy[x].direction[i].blue;
1609       if (image->colorspace == CMYKColorspace)
1610         channel_features[BlackPixelChannel].sum_variance[i]+=
1611           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1612           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1613           density_xy[x].direction[i].black;
1614       if (image->alpha_trait == BlendPixelTrait)
1615         channel_features[AlphaPixelChannel].sum_variance[i]+=
1616           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1617           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1618           density_xy[x].direction[i].alpha;
1619     }
1620   }
1621   /*
1622     Compute more texture features.
1623   */
1624 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1625   #pragma omp parallel for schedule(static,4) shared(status) \
1626     magick_threads(image,image,number_grays,1)
1627 #endif
1628   for (i=0; i < 4; i++)
1629   {
1630     register ssize_t
1631       y;
1632
1633     for (y=0; y < (ssize_t) number_grays; y++)
1634     {
1635       register ssize_t
1636         x;
1637
1638       for (x=0; x < (ssize_t) number_grays; x++)
1639       {
1640         /*
1641           Sum of Squares: Variance
1642         */
1643         variance.direction[i].red+=(y-mean.direction[i].red+1)*
1644           (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1645         variance.direction[i].green+=(y-mean.direction[i].green+1)*
1646           (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1647         variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1648           (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1649         if (image->colorspace == CMYKColorspace)
1650           variance.direction[i].black+=(y-mean.direction[i].black+1)*
1651             (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1652         if (image->alpha_trait == BlendPixelTrait)
1653           variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1654             (y-mean.direction[i].alpha+1)*
1655             cooccurrence[x][y].direction[i].alpha;
1656         /*
1657           Sum average / Difference Variance.
1658         */
1659         density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1660           cooccurrence[x][y].direction[i].red;
1661         density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1662           cooccurrence[x][y].direction[i].green;
1663         density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1664           cooccurrence[x][y].direction[i].blue;
1665         if (image->colorspace == CMYKColorspace)
1666           density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1667             cooccurrence[x][y].direction[i].black;
1668         if (image->alpha_trait == BlendPixelTrait)
1669           density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1670             cooccurrence[x][y].direction[i].alpha;
1671         /*
1672           Information Measures of Correlation.
1673         */
1674         entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1675           MagickLog10(cooccurrence[x][y].direction[i].red);
1676         entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1677           MagickLog10(cooccurrence[x][y].direction[i].green);
1678         entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1679           MagickLog10(cooccurrence[x][y].direction[i].blue);
1680         if (image->colorspace == CMYKColorspace)
1681           entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1682             MagickLog10(cooccurrence[x][y].direction[i].black);
1683         if (image->alpha_trait == BlendPixelTrait)
1684           entropy_xy.direction[i].alpha-=
1685             cooccurrence[x][y].direction[i].alpha*MagickLog10(
1686             cooccurrence[x][y].direction[i].alpha);
1687         entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1688           MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
1689         entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1690           MagickLog10(density_x[x].direction[i].green*
1691           density_y[y].direction[i].green));
1692         entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1693           MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1694         if (image->colorspace == CMYKColorspace)
1695           entropy_xy1.direction[i].black-=(
1696             cooccurrence[x][y].direction[i].black*MagickLog10(
1697             density_x[x].direction[i].black*density_y[y].direction[i].black));
1698         if (image->alpha_trait == BlendPixelTrait)
1699           entropy_xy1.direction[i].alpha-=(
1700             cooccurrence[x][y].direction[i].alpha*MagickLog10(
1701             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1702         entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1703           density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
1704           density_y[y].direction[i].red));
1705         entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1706           density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
1707           density_y[y].direction[i].green));
1708         entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1709           density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
1710           density_y[y].direction[i].blue));
1711         if (image->colorspace == CMYKColorspace)
1712           entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1713             density_y[y].direction[i].black*MagickLog10(
1714             density_x[x].direction[i].black*density_y[y].direction[i].black));
1715         if (image->alpha_trait == BlendPixelTrait)
1716           entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1717             density_y[y].direction[i].alpha*MagickLog10(
1718             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1719       }
1720     }
1721     channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1722       variance.direction[i].red;
1723     channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1724       variance.direction[i].green;
1725     channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1726       variance.direction[i].blue;
1727     if (image->colorspace == CMYKColorspace)
1728       channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1729         variance.direction[i].black;
1730     if (image->alpha_trait == BlendPixelTrait)
1731       channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1732         variance.direction[i].alpha;
1733   }
1734   /*
1735     Compute more texture features.
1736   */
1737   (void) ResetMagickMemory(&variance,0,sizeof(variance));
1738   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
1739 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1740   #pragma omp parallel for schedule(static,4) shared(status) \
1741     magick_threads(image,image,number_grays,1)
1742 #endif
1743   for (i=0; i < 4; i++)
1744   {
1745     register ssize_t
1746       x;
1747
1748     for (x=0; x < (ssize_t) number_grays; x++)
1749     {
1750       /*
1751         Difference variance.
1752       */
1753       variance.direction[i].red+=density_xy[x].direction[i].red;
1754       variance.direction[i].green+=density_xy[x].direction[i].green;
1755       variance.direction[i].blue+=density_xy[x].direction[i].blue;
1756       if (image->colorspace == CMYKColorspace)
1757         variance.direction[i].black+=density_xy[x].direction[i].black;
1758       if (image->alpha_trait == BlendPixelTrait)
1759         variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1760       sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1761         density_xy[x].direction[i].red;
1762       sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1763         density_xy[x].direction[i].green;
1764       sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1765         density_xy[x].direction[i].blue;
1766       if (image->colorspace == CMYKColorspace)
1767         sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1768           density_xy[x].direction[i].black;
1769       if (image->alpha_trait == BlendPixelTrait)
1770         sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1771           density_xy[x].direction[i].alpha;
1772       /*
1773         Difference entropy.
1774       */
1775       channel_features[RedPixelChannel].difference_entropy[i]-=
1776         density_xy[x].direction[i].red*
1777         MagickLog10(density_xy[x].direction[i].red);
1778       channel_features[GreenPixelChannel].difference_entropy[i]-=
1779         density_xy[x].direction[i].green*
1780         MagickLog10(density_xy[x].direction[i].green);
1781       channel_features[BluePixelChannel].difference_entropy[i]-=
1782         density_xy[x].direction[i].blue*
1783         MagickLog10(density_xy[x].direction[i].blue);
1784       if (image->colorspace == CMYKColorspace)
1785         channel_features[BlackPixelChannel].difference_entropy[i]-=
1786           density_xy[x].direction[i].black*
1787           MagickLog10(density_xy[x].direction[i].black);
1788       if (image->alpha_trait == BlendPixelTrait)
1789         channel_features[AlphaPixelChannel].difference_entropy[i]-=
1790           density_xy[x].direction[i].alpha*
1791           MagickLog10(density_xy[x].direction[i].alpha);
1792       /*
1793         Information Measures of Correlation.
1794       */
1795       entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1796         MagickLog10(density_x[x].direction[i].red));
1797       entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1798         MagickLog10(density_x[x].direction[i].green));
1799       entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1800         MagickLog10(density_x[x].direction[i].blue));
1801       if (image->colorspace == CMYKColorspace)
1802         entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1803           MagickLog10(density_x[x].direction[i].black));
1804       if (image->alpha_trait == BlendPixelTrait)
1805         entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1806           MagickLog10(density_x[x].direction[i].alpha));
1807       entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1808         MagickLog10(density_y[x].direction[i].red));
1809       entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1810         MagickLog10(density_y[x].direction[i].green));
1811       entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1812         MagickLog10(density_y[x].direction[i].blue));
1813       if (image->colorspace == CMYKColorspace)
1814         entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1815           MagickLog10(density_y[x].direction[i].black));
1816       if (image->alpha_trait == BlendPixelTrait)
1817         entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1818           MagickLog10(density_y[x].direction[i].alpha));
1819     }
1820     /*
1821       Difference variance.
1822     */
1823     channel_features[RedPixelChannel].difference_variance[i]=
1824       (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1825       (variance.direction[i].red*variance.direction[i].red))/
1826       ((double) number_grays*number_grays*number_grays*number_grays);
1827     channel_features[GreenPixelChannel].difference_variance[i]=
1828       (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1829       (variance.direction[i].green*variance.direction[i].green))/
1830       ((double) number_grays*number_grays*number_grays*number_grays);
1831     channel_features[BluePixelChannel].difference_variance[i]=
1832       (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1833       (variance.direction[i].blue*variance.direction[i].blue))/
1834       ((double) number_grays*number_grays*number_grays*number_grays);
1835     if (image->colorspace == CMYKColorspace)
1836       channel_features[BlackPixelChannel].difference_variance[i]=
1837         (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1838         (variance.direction[i].black*variance.direction[i].black))/
1839         ((double) number_grays*number_grays*number_grays*number_grays);
1840     if (image->alpha_trait == BlendPixelTrait)
1841       channel_features[AlphaPixelChannel].difference_variance[i]=
1842         (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1843         (variance.direction[i].alpha*variance.direction[i].alpha))/
1844         ((double) number_grays*number_grays*number_grays*number_grays);
1845     /*
1846       Information Measures of Correlation.
1847     */
1848     channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1849       (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1850       (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1851        entropy_x.direction[i].red : entropy_y.direction[i].red);
1852     channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1853       (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1854       (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1855        entropy_x.direction[i].green : entropy_y.direction[i].green);
1856     channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1857       (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1858       (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1859        entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1860     if (image->colorspace == CMYKColorspace)
1861       channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1862         (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1863         (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1864          entropy_x.direction[i].black : entropy_y.direction[i].black);
1865     if (image->alpha_trait == BlendPixelTrait)
1866       channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1867         (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1868         (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1869          entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1870     channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1871       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
1872       entropy_xy.direction[i].red)))));
1873     channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1874       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
1875       entropy_xy.direction[i].green)))));
1876     channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1877       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
1878       entropy_xy.direction[i].blue)))));
1879     if (image->colorspace == CMYKColorspace)
1880       channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1881         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black-
1882         entropy_xy.direction[i].black)))));
1883     if (image->alpha_trait == BlendPixelTrait)
1884       channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1885         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha-
1886         entropy_xy.direction[i].alpha)))));
1887   }
1888   /*
1889     Compute more texture features.
1890   */
1891 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1892   #pragma omp parallel for schedule(static,4) shared(status) \
1893     magick_threads(image,image,number_grays,1)
1894 #endif
1895   for (i=0; i < 4; i++)
1896   {
1897     ssize_t
1898       z;
1899
1900     for (z=0; z < (ssize_t) number_grays; z++)
1901     {
1902       register ssize_t
1903         y;
1904
1905       ChannelStatistics
1906         pixel;
1907
1908       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
1909       for (y=0; y < (ssize_t) number_grays; y++)
1910       {
1911         register ssize_t
1912           x;
1913
1914         for (x=0; x < (ssize_t) number_grays; x++)
1915         {
1916           /*
1917             Contrast:  amount of local variations present in an image.
1918           */
1919           if (((y-x) == z) || ((x-y) == z))
1920             {
1921               pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1922               pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1923               pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1924               if (image->colorspace == CMYKColorspace)
1925                 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1926               if (image->alpha_trait == BlendPixelTrait)
1927                 pixel.direction[i].alpha+=
1928                   cooccurrence[x][y].direction[i].alpha;
1929             }
1930           /*
1931             Maximum Correlation Coefficient.
1932           */
1933           Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1934             cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1935             density_y[x].direction[i].red;
1936           Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1937             cooccurrence[y][x].direction[i].green/
1938             density_x[z].direction[i].green/density_y[x].direction[i].red;
1939           Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1940             cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
1941             density_y[x].direction[i].blue;
1942           if (image->colorspace == CMYKColorspace)
1943             Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1944               cooccurrence[y][x].direction[i].black/
1945               density_x[z].direction[i].black/density_y[x].direction[i].black;
1946           if (image->alpha_trait == BlendPixelTrait)
1947             Q[z][y].direction[i].alpha+=
1948               cooccurrence[z][x].direction[i].alpha*
1949               cooccurrence[y][x].direction[i].alpha/
1950               density_x[z].direction[i].alpha/
1951               density_y[x].direction[i].alpha;
1952         }
1953       }
1954       channel_features[RedPixelChannel].contrast[i]+=z*z*
1955         pixel.direction[i].red;
1956       channel_features[GreenPixelChannel].contrast[i]+=z*z*
1957         pixel.direction[i].green;
1958       channel_features[BluePixelChannel].contrast[i]+=z*z*
1959         pixel.direction[i].blue;
1960       if (image->colorspace == CMYKColorspace)
1961         channel_features[BlackPixelChannel].contrast[i]+=z*z*
1962           pixel.direction[i].black;
1963       if (image->alpha_trait == BlendPixelTrait)
1964         channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1965           pixel.direction[i].alpha;
1966     }
1967     /*
1968       Maximum Correlation Coefficient.
1969       Future: return second largest eigenvalue of Q.
1970     */
1971     channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1972       sqrt((double) -1.0);
1973     channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1974       sqrt((double) -1.0);
1975     channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1976       sqrt((double) -1.0);
1977     if (image->colorspace == CMYKColorspace)
1978       channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1979         sqrt((double) -1.0);
1980     if (image->alpha_trait == BlendPixelTrait)
1981       channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1982         sqrt((double) -1.0);
1983   }
1984   /*
1985     Relinquish resources.
1986   */
1987   sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1988   for (i=0; i < (ssize_t) number_grays; i++)
1989     Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1990   Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1991   density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1992   density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1993   density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1994   for (i=0; i < (ssize_t) number_grays; i++)
1995     cooccurrence[i]=(ChannelStatistics *)
1996       RelinquishMagickMemory(cooccurrence[i]);
1997   cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1998   return(channel_features);
1999 }