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

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %               SSSSS  EEEEE   GGGG  M   M  EEEEE  N   N  TTTTT               %
   7 %               SS     E      G      MM MM  E      NN  N    T                 %
   8 %                SSS   EEE    G GGG  M M M  EEE    N N N    T                 %
   9 %                  SS  E      G   G  M   M  E      N  NN    T                 %
  10 %               SSSSS  EEEEE   GGGG  M   M  EEEEE  N   N    T                 %
  11 %                                                                             %
  12 %                                                                             %
  13 %    MagickCore Methods to Segment an Image with Thresholding Fuzzy c-Means   %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                                John Cristy                                  %
  17 %                                April 1993                                   %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2011 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 %  Segment segments an image by analyzing the histograms of the color
  37 %  components and identifying units that are homogeneous with the fuzzy
  38 %  c-means technique.  The scale-space filter analyzes the histograms of
  39 %  the three color components of the image and identifies a set of
  40 %  classes.  The extents of each class is used to coarsely segment the
  41 %  image with thresholding.  The color associated with each class is
  42 %  determined by the mean color of all pixels within the extents of a
  43 %  particular class.  Finally, any unclassified pixels are assigned to
  44 %  the closest class with the fuzzy c-means technique.
  45 %
  46 %  The fuzzy c-Means algorithm can be summarized as follows:
  47 %
  48 %    o Build a histogram, one for each color component of the image.
  49 %
  50 %    o For each histogram, successively apply the scale-space filter and
  51 %      build an interval tree of zero crossings in the second derivative
  52 %      at each scale.  Analyze this scale-space ``fingerprint'' to
  53 %      determine which peaks and valleys in the histogram are most
  54 %      predominant.
  55 %
  56 %    o The fingerprint defines intervals on the axis of the histogram.
  57 %      Each interval contains either a minima or a maxima in the original
  58 %      signal.  If each color component lies within the maxima interval,
  59 %      that pixel is considered ``classified'' and is assigned an unique
  60 %      class number.
  61 %
  62 %    o Any pixel that fails to be classified in the above thresholding
  63 %      pass is classified using the fuzzy c-Means technique.  It is
  64 %      assigned to one of the classes discovered in the histogram analysis
  65 %      phase.
  66 %
  67 %  The fuzzy c-Means technique attempts to cluster a pixel by finding
  68 %  the local minima of the generalized within group sum of squared error
  69 %  objective function.  A pixel is assigned to the closest class of
  70 %  which the fuzzy membership has a maximum value.
  71 %
  72 %  Segment is strongly based on software written by Andy Gallo,
  73 %  University of Delaware.
  74 %
  75 %  The following reference was used in creating this program:
  76 %
  77 %    Young Won Lim, Sang Uk Lee, "On The Color Image Segmentation
  78 %    Algorithm Based on the Thresholding and the Fuzzy c-Means
  79 %    Techniques", Pattern Recognition, Volume 23, Number 9, pages
  80 %    935-952, 1990.
  81 %
  82 %
  83 */
  84 \f
  85 #include "MagickCore/studio.h"
  86 #include "MagickCore/cache.h"
  87 #include "MagickCore/color.h"
  88 #include "MagickCore/colormap.h"
  89 #include "MagickCore/colorspace.h"
  90 #include "MagickCore/colorspace-private.h"
  91 #include "MagickCore/exception.h"
  92 #include "MagickCore/exception-private.h"
  93 #include "MagickCore/image.h"
  94 #include "MagickCore/image-private.h"
  95 #include "MagickCore/memory_.h"
  96 #include "MagickCore/monitor.h"
  97 #include "MagickCore/monitor-private.h"
  98 #include "MagickCore/pixel-accessor.h"
  99 #include "MagickCore/quantize.h"
 100 #include "MagickCore/quantum.h"
 101 #include "MagickCore/quantum-private.h"
 102 #include "MagickCore/segment.h"
 103 #include "MagickCore/string_.h"
 104 \f
 105 /*
 106   Define declarations.
 107 */
 108 #define MaxDimension  3
 109 #define DeltaTau  0.5f
 110 #if defined(FastClassify)
 111 #define WeightingExponent  2.0
 112 #define SegmentPower(ratio) (ratio)
 113 #else
 114 #define WeightingExponent  2.5
 115 #define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0)));
 116 #endif
 117 #define Tau  5.2f
 118 \f
 119 /*
 120   Typedef declarations.
 121 */
 122 typedef struct _ExtentPacket
 123 {
 124   MagickRealType
 125     center;
 126
 127   ssize_t
 128     index,
 129     left,
 130     right;
 131 } ExtentPacket;
 132
 133 typedef struct _Cluster
 134 {
 135   struct _Cluster
 136     *next;
 137
 138   ExtentPacket
 139     red,
 140     green,
 141     blue;
 142
 143   ssize_t
 144     count,
 145     id;
 146 } Cluster;
 147
 148 typedef struct _IntervalTree
 149 {
 150   MagickRealType
 151     tau;
 152
 153   ssize_t
 154     left,
 155     right;
 156
 157   MagickRealType
 158     mean_stability,
 159     stability;
 160
 161   struct _IntervalTree
 162     *sibling,
 163     *child;
 164 } IntervalTree;
 165
 166 typedef struct _ZeroCrossing
 167 {
 168   MagickRealType
 169     tau,
 170     histogram[256];
 171
 172   short
 173     crossings[256];
 174 } ZeroCrossing;
 175 \f
 176 /*
 177   Constant declarations.
 178 */
 179 static const int
 180   Blue = 2,
 181   Green = 1,
 182   Red = 0,
 183   SafeMargin = 3,
 184   TreeLength = 600;
 185 \f
 186 /*
 187   Method prototypes.
 188 */
 189 static MagickRealType
 190   OptimalTau(const ssize_t *,const double,const double,const double,
 191     const double,short *);
 192
 193 static ssize_t
 194   DefineRegion(const short *,ExtentPacket *);
 195
 196 static void
 197   InitializeHistogram(const Image *,ssize_t **,ExceptionInfo *),
 198   ScaleSpace(const ssize_t *,const MagickRealType,MagickRealType *),
 199   ZeroCrossHistogram(MagickRealType *,const MagickRealType,short *);
 200 \f
 201 /*
 202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 203 %                                                                             %
 204 %                                                                             %
 205 %                                                                             %
 206 +   C l a s s i f y                                                           %
 207 %                                                                             %
 208 %                                                                             %
 209 %                                                                             %
 210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 211 %
 212 %  Classify() defines one or more classes.  Each pixel is thresholded to
 213 %  determine which class it belongs to.  If the class is not identified it is
 214 %  assigned to the closest class based on the fuzzy c-Means technique.
 215 %
 216 %  The format of the Classify method is:
 217 %
 218 %      MagickBooleanType Classify(Image *image,short **extrema,
 219 %        const MagickRealType cluster_threshold,
 220 %        const MagickRealType weighting_exponent,
 221 %        const MagickBooleanType verbose)
 222 %
 223 %  A description of each parameter follows.
 224 %
 225 %    o image: the image.
 226 %
 227 %    o extrema:  Specifies a pointer to an array of integers.  They
 228 %      represent the peaks and valleys of the histogram for each color
 229 %      component.
 230 %
 231 %    o cluster_threshold:  This MagickRealType represents the minimum number of
 232 %      pixels contained in a hexahedra before it can be considered valid
 233 %      (expressed as a percentage).
 234 %
 235 %    o weighting_exponent: Specifies the membership weighting exponent.
 236 %
 237 %    o verbose:  A value greater than zero prints detailed information about
 238 %      the identified classes.
 239 %
 240 */
 241 static MagickBooleanType Classify(Image *image,short **extrema,
 242   const MagickRealType cluster_threshold,
 243   const MagickRealType weighting_exponent,const MagickBooleanType verbose)
 244 {
 245 #define SegmentImageTag  "Segment/Image"
 246
 247   CacheView
 248     *image_view;
 249
 250   Cluster
 251     *cluster,
 252     *head,
 253     *last_cluster,
 254     *next_cluster;
 255
 256   ExceptionInfo
 257     *exception;
 258
 259   ExtentPacket
 260     blue,
 261     green,
 262     red;
 263
 264   MagickOffsetType
 265     progress;
 266
 267   MagickRealType
 268     *free_squares;
 269
 270   MagickStatusType
 271     status;
 272
 273   register ssize_t
 274     i;
 275
 276   register MagickRealType
 277     *squares;
 278
 279   size_t
 280     number_clusters;
 281
 282   ssize_t
 283     count,
 284     y;
 285
 286   /*
 287     Form clusters.
 288   */
 289   cluster=(Cluster *) NULL;
 290   head=(Cluster *) NULL;
 291   (void) ResetMagickMemory(&red,0,sizeof(red));
 292   (void) ResetMagickMemory(&green,0,sizeof(green));
 293   (void) ResetMagickMemory(&blue,0,sizeof(blue));
 294   while (DefineRegion(extrema[Red],&red) != 0)
 295   {
 296     green.index=0;
 297     while (DefineRegion(extrema[Green],&green) != 0)
 298     {
 299       blue.index=0;
 300       while (DefineRegion(extrema[Blue],&blue) != 0)
 301       {
 302         /*
 303           Allocate a new class.
 304         */
 305         if (head != (Cluster *) NULL)
 306           {
 307             cluster->next=(Cluster *) AcquireMagickMemory(
 308               sizeof(*cluster->next));
 309             cluster=cluster->next;
 310           }
 311         else
 312           {
 313             cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
 314             head=cluster;
 315           }
 316         if (cluster == (Cluster *) NULL)
 317           ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
 318             image->filename);
 319         /*
 320           Initialize a new class.
 321         */
 322         cluster->count=0;
 323         cluster->red=red;
 324         cluster->green=green;
 325         cluster->blue=blue;
 326         cluster->next=(Cluster *) NULL;
 327       }
 328     }
 329   }
 330   if (head == (Cluster *) NULL)
 331     {
 332       /*
 333         No classes were identified-- create one.
 334       */
 335       cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
 336       if (cluster == (Cluster *) NULL)
 337         ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
 338           image->filename);
 339       /*
 340         Initialize a new class.
 341       */
 342       cluster->count=0;
 343       cluster->red=red;
 344       cluster->green=green;
 345       cluster->blue=blue;
 346       cluster->next=(Cluster *) NULL;
 347       head=cluster;
 348     }
 349   /*
 350     Count the pixels for each cluster.
 351   */
 352   status=MagickTrue;
 353   count=0;
 354   progress=0;
 355   exception=(&image->exception);
 356   image_view=AcquireCacheView(image);
 357   for (y=0; y < (ssize_t) image->rows; y++)
 358   {
 359     register const Quantum
 360       *p;
 361
 362     register ssize_t
 363       x;
 364
 365     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 366     if (p == (const Quantum *) NULL)
 367       break;
 368     for (x=0; x < (ssize_t) image->columns; x++)
 369     {
 370       for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
 371         if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >=
 372              (cluster->red.left-SafeMargin)) &&
 373             ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <=
 374              (cluster->red.right+SafeMargin)) &&
 375             ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >=
 376              (cluster->green.left-SafeMargin)) &&
 377             ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <=
 378              (cluster->green.right+SafeMargin)) &&
 379             ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >=
 380              (cluster->blue.left-SafeMargin)) &&
 381             ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <=
 382              (cluster->blue.right+SafeMargin)))
 383           {
 384             /*
 385               Count this pixel.
 386             */
 387             count++;
 388             cluster->red.center+=(MagickRealType) ScaleQuantumToChar(
 389               GetPixelRed(image,p));
 390             cluster->green.center+=(MagickRealType) ScaleQuantumToChar(
 391               GetPixelGreen(image,p));
 392             cluster->blue.center+=(MagickRealType) ScaleQuantumToChar(
 393               GetPixelBlue(image,p));
 394             cluster->count++;
 395             break;
 396           }
 397       p+=GetPixelChannels(image);
 398     }
 399     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 400       {
 401         MagickBooleanType
 402           proceed;
 403
 404 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 405         #pragma omp critical (MagickCore_Classify)
 406 #endif
 407         proceed=SetImageProgress(image,SegmentImageTag,progress++,
 408           2*image->rows);
 409         if (proceed == MagickFalse)
 410           status=MagickFalse;
 411       }
 412   }
 413   image_view=DestroyCacheView(image_view);
 414   /*
 415     Remove clusters that do not meet minimum cluster threshold.
 416   */
 417   count=0;
 418   last_cluster=head;
 419   next_cluster=head;
 420   for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
 421   {
 422     next_cluster=cluster->next;
 423     if ((cluster->count > 0) &&
 424         (cluster->count >= (count*cluster_threshold/100.0)))
 425       {
 426         /*
 427           Initialize cluster.
 428         */
 429         cluster->id=count;
 430         cluster->red.center/=cluster->count;
 431         cluster->green.center/=cluster->count;
 432         cluster->blue.center/=cluster->count;
 433         count++;
 434         last_cluster=cluster;
 435         continue;
 436       }
 437     /*
 438       Delete cluster.
 439     */
 440     if (cluster == head)
 441       head=next_cluster;
 442     else
 443       last_cluster->next=next_cluster;
 444     cluster=(Cluster *) RelinquishMagickMemory(cluster);
 445   }
 446   number_clusters=(size_t) count;
 447   if (verbose != MagickFalse)
 448     {
 449       /*
 450         Print cluster statistics.
 451       */
 452       (void) FormatLocaleFile(stdout,"Fuzzy C-means Statistics\n");
 453       (void) FormatLocaleFile(stdout,"===================\n\n");
 454       (void) FormatLocaleFile(stdout,"\tCluster Threshold = %g\n",(double)
 455         cluster_threshold);
 456       (void) FormatLocaleFile(stdout,"\tWeighting Exponent = %g\n",(double)
 457         weighting_exponent);
 458       (void) FormatLocaleFile(stdout,"\tTotal Number of Clusters = %.20g\n\n",
 459         (double) number_clusters);
 460       /*
 461         Print the total number of points per cluster.
 462       */
 463       (void) FormatLocaleFile(stdout,"\n\nNumber of Vectors Per Cluster\n");
 464       (void) FormatLocaleFile(stdout,"=============================\n\n");
 465       for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
 466         (void) FormatLocaleFile(stdout,"Cluster #%.20g = %.20g\n",(double)
 467           cluster->id,(double) cluster->count);
 468       /*
 469         Print the cluster extents.
 470       */
 471       (void) FormatLocaleFile(stdout,
 472         "\n\n\nCluster Extents:        (Vector Size: %d)\n",MaxDimension);
 473       (void) FormatLocaleFile(stdout,"================");
 474       for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
 475       {
 476         (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
 477           cluster->id);
 478         (void) FormatLocaleFile(stdout,
 479           "%.20g-%.20g  %.20g-%.20g  %.20g-%.20g\n",(double)
 480           cluster->red.left,(double) cluster->red.right,(double)
 481           cluster->green.left,(double) cluster->green.right,(double)
 482           cluster->blue.left,(double) cluster->blue.right);
 483       }
 484       /*
 485         Print the cluster center values.
 486       */
 487       (void) FormatLocaleFile(stdout,
 488         "\n\n\nCluster Center Values:        (Vector Size: %d)\n",MaxDimension);
 489       (void) FormatLocaleFile(stdout,"=====================");
 490       for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
 491       {
 492         (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
 493           cluster->id);
 494         (void) FormatLocaleFile(stdout,"%g  %g  %g\n",(double)
 495           cluster->red.center,(double) cluster->green.center,(double)
 496           cluster->blue.center);
 497       }
 498       (void) FormatLocaleFile(stdout,"\n");
 499     }
 500   if (number_clusters > 256)
 501     ThrowBinaryException(ImageError,"TooManyClusters",image->filename);
 502   /*
 503     Speed up distance calculations.
 504   */
 505   squares=(MagickRealType *) AcquireQuantumMemory(513UL,sizeof(*squares));
 506   if (squares == (MagickRealType *) NULL)
 507     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
 508       image->filename);
 509   squares+=255;
 510   for (i=(-255); i <= 255; i++)
 511     squares[i]=(MagickRealType) i*(MagickRealType) i;
 512   /*
 513     Allocate image colormap.
 514   */
 515   if (AcquireImageColormap(image,number_clusters) == MagickFalse)
 516     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
 517       image->filename);
 518   i=0;
 519   for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
 520   {
 521     image->colormap[i].red=ScaleCharToQuantum((unsigned char)
 522       (cluster->red.center+0.5));
 523     image->colormap[i].green=ScaleCharToQuantum((unsigned char)
 524       (cluster->green.center+0.5));
 525     image->colormap[i].blue=ScaleCharToQuantum((unsigned char)
 526       (cluster->blue.center+0.5));
 527     i++;
 528   }
 529   /*
 530     Do course grain classes.
 531   */
 532   exception=(&image->exception);
 533   image_view=AcquireCacheView(image);
 534 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 535   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 536 #endif
 537   for (y=0; y < (ssize_t) image->rows; y++)
 538   {
 539     Cluster
 540       *cluster;
 541
 542     register const PixelPacket
 543       *restrict p;
 544
 545     register ssize_t
 546       x;
 547
 548     register Quantum
 549       *restrict q;
 550
 551     if (status == MagickFalse)
 552       continue;
 553     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 554     if (q == (const Quantum *) NULL)
 555       {
 556         status=MagickFalse;
 557         continue;
 558       }
 559     for (x=0; x < (ssize_t) image->columns; x++)
 560     {
 561       SetPixelIndex(image,0,q);
 562       for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
 563       {
 564         if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,q)) >=
 565              (cluster->red.left-SafeMargin)) &&
 566             ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,q)) <=
 567              (cluster->red.right+SafeMargin)) &&
 568             ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q)) >=
 569              (cluster->green.left-SafeMargin)) &&
 570             ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q)) <=
 571              (cluster->green.right+SafeMargin)) &&
 572             ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q)) >=
 573              (cluster->blue.left-SafeMargin)) &&
 574             ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q)) <=
 575              (cluster->blue.right+SafeMargin)))
 576           {
 577             /*
 578               Classify this pixel.
 579             */
 580             SetPixelIndex(image,(Quantum) cluster->id,q);
 581             break;
 582           }
 583       }
 584       if (cluster == (Cluster *) NULL)
 585         {
 586           MagickRealType
 587             distance_squared,
 588             local_minima,
 589             numerator,
 590             ratio,
 591             sum;
 592
 593           register ssize_t
 594             j,
 595             k;
 596
 597           /*
 598             Compute fuzzy membership.
 599           */
 600           local_minima=0.0;
 601           for (j=0; j < (ssize_t) image->colors; j++)
 602           {
 603             sum=0.0;
 604             p=image->colormap+j;
 605             distance_squared=squares[(ssize_t) ScaleQuantumToChar(
 606               GetPixelRed(image,q))-(ssize_t)
 607               ScaleQuantumToChar(p->red)]+squares[(ssize_t)
 608               ScaleQuantumToChar(GetPixelGreen(image,q))-(ssize_t)
 609               ScaleQuantumToChar(p->green)]+squares[(ssize_t)
 610               ScaleQuantumToChar(GetPixelBlue(image,q))-(ssize_t)
 611               ScaleQuantumToChar(p->blue)];
 612             numerator=distance_squared;
 613             for (k=0; k < (ssize_t) image->colors; k++)
 614             {
 615               p=image->colormap+k;
 616                 distance_squared=squares[(ssize_t) ScaleQuantumToChar(
 617                   GetPixelRed(image,q))-(ssize_t)
 618                   ScaleQuantumToChar(p->red)]+squares[(ssize_t)
 619                   ScaleQuantumToChar(GetPixelGreen(image,q))-(ssize_t)
 620                   ScaleQuantumToChar(p->green)]+squares[(ssize_t)
 621                   ScaleQuantumToChar(GetPixelBlue(image,q))-(ssize_t)
 622                   ScaleQuantumToChar(p->blue)];
 623               ratio=numerator/distance_squared;
 624               sum+=SegmentPower(ratio);
 625             }
 626             if ((sum != 0.0) && ((1.0/sum) > local_minima))
 627               {
 628                 /*
 629                   Classify this pixel.
 630                 */
 631                 local_minima=1.0/sum;
 632                 SetPixelIndex(image,(Quantum) j,q);
 633               }
 634           }
 635         }
 636       q+=GetPixelChannels(image);
 637     }
 638     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 639       status=MagickFalse;
 640     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 641       {
 642         MagickBooleanType
 643           proceed;
 644
 645 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 646         #pragma omp critical (MagickCore_Classify)
 647 #endif
 648         proceed=SetImageProgress(image,SegmentImageTag,progress++,
 649           2*image->rows);
 650         if (proceed == MagickFalse)
 651           status=MagickFalse;
 652       }
 653   }
 654   image_view=DestroyCacheView(image_view);
 655   status&=SyncImage(image);
 656   /*
 657     Relinquish resources.
 658   */
 659   for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
 660   {
 661     next_cluster=cluster->next;
 662     cluster=(Cluster *) RelinquishMagickMemory(cluster);
 663   }
 664   squares-=255;
 665   free_squares=squares;
 666   free_squares=(MagickRealType *) RelinquishMagickMemory(free_squares);
 667   return(MagickTrue);
 668 }
 669 \f
 670 /*
 671 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 672 %                                                                             %
 673 %                                                                             %
 674 %                                                                             %
 675 +   C o n s o l i d a t e C r o s s i n g s                                   %
 676 %                                                                             %
 677 %                                                                             %
 678 %                                                                             %
 679 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 680 %
 681 %  ConsolidateCrossings() guarantees that an even number of zero crossings
 682 %  always lie between two crossings.
 683 %
 684 %  The format of the ConsolidateCrossings method is:
 685 %
 686 %      ConsolidateCrossings(ZeroCrossing *zero_crossing,
 687 %        const size_t number_crossings)
 688 %
 689 %  A description of each parameter follows.
 690 %
 691 %    o zero_crossing: Specifies an array of structures of type ZeroCrossing.
 692 %
 693 %    o number_crossings: This size_t specifies the number of elements
 694 %      in the zero_crossing array.
 695 %
 696 */
 697
 698 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
 699 {
 700   if (x < 0)
 701     return(-x);
 702   return(x);
 703 }
 704
 705 static inline ssize_t MagickMax(const ssize_t x,const ssize_t y)
 706 {
 707   if (x > y)
 708     return(x);
 709   return(y);
 710 }
 711
 712 static inline ssize_t MagickMin(const ssize_t x,const ssize_t y)
 713 {
 714   if (x < y)
 715     return(x);
 716   return(y);
 717 }
 718
 719 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
 720   const size_t number_crossings)
 721 {
 722   register ssize_t
 723     i,
 724     j,
 725     k,
 726     l;
 727
 728   ssize_t
 729     center,
 730     correct,
 731     count,
 732     left,
 733     right;
 734
 735   /*
 736     Consolidate zero crossings.
 737   */
 738   for (i=(ssize_t) number_crossings-1; i >= 0; i--)
 739     for (j=0; j <= 255; j++)
 740     {
 741       if (zero_crossing[i].crossings[j] == 0)
 742         continue;
 743       /*
 744         Find the entry that is closest to j and still preserves the
 745         property that there are an even number of crossings between
 746         intervals.
 747       */
 748       for (k=j-1; k > 0; k--)
 749         if (zero_crossing[i+1].crossings[k] != 0)
 750           break;
 751       left=MagickMax(k,0);
 752       center=j;
 753       for (k=j+1; k < 255; k++)
 754         if (zero_crossing[i+1].crossings[k] != 0)
 755           break;
 756       right=MagickMin(k,255);
 757       /*
 758         K is the zero crossing just left of j.
 759       */
 760       for (k=j-1; k > 0; k--)
 761         if (zero_crossing[i].crossings[k] != 0)
 762           break;
 763       if (k < 0)
 764         k=0;
 765       /*
 766         Check center for an even number of crossings between k and j.
 767       */
 768       correct=(-1);
 769       if (zero_crossing[i+1].crossings[j] != 0)
 770         {
 771           count=0;
 772           for (l=k+1; l < center; l++)
 773             if (zero_crossing[i+1].crossings[l] != 0)
 774               count++;
 775           if (((count % 2) == 0) && (center != k))
 776             correct=center;
 777         }
 778       /*
 779         Check left for an even number of crossings between k and j.
 780       */
 781       if (correct == -1)
 782         {
 783           count=0;
 784           for (l=k+1; l < left; l++)
 785             if (zero_crossing[i+1].crossings[l] != 0)
 786               count++;
 787           if (((count % 2) == 0) && (left != k))
 788             correct=left;
 789         }
 790       /*
 791         Check right for an even number of crossings between k and j.
 792       */
 793       if (correct == -1)
 794         {
 795           count=0;
 796           for (l=k+1; l < right; l++)
 797             if (zero_crossing[i+1].crossings[l] != 0)
 798               count++;
 799           if (((count % 2) == 0) && (right != k))
 800             correct=right;
 801         }
 802       l=(ssize_t) zero_crossing[i].crossings[j];
 803       zero_crossing[i].crossings[j]=0;
 804       if (correct != -1)
 805         zero_crossing[i].crossings[correct]=(short) l;
 806     }
 807 }
 808 \f
 809 /*
 810 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 811 %                                                                             %
 812 %                                                                             %
 813 %                                                                             %
 814 +   D e f i n e R e g i o n                                                   %
 815 %                                                                             %
 816 %                                                                             %
 817 %                                                                             %
 818 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 819 %
 820 %  DefineRegion() defines the left and right boundaries of a peak region.
 821 %
 822 %  The format of the DefineRegion method is:
 823 %
 824 %      ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
 825 %
 826 %  A description of each parameter follows.
 827 %
 828 %    o extrema:  Specifies a pointer to an array of integers.  They
 829 %      represent the peaks and valleys of the histogram for each color
 830 %      component.
 831 %
 832 %    o extents:  This pointer to an ExtentPacket represent the extends
 833 %      of a particular peak or valley of a color component.
 834 %
 835 */
 836 static ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
 837 {
 838   /*
 839     Initialize to default values.
 840   */
 841   extents->left=0;
 842   extents->center=0.0;
 843   extents->right=255;
 844   /*
 845     Find the left side (maxima).
 846   */
 847   for ( ; extents->index <= 255; extents->index++)
 848     if (extrema[extents->index] > 0)
 849       break;
 850   if (extents->index > 255)
 851     return(MagickFalse);  /* no left side - no region exists */
 852   extents->left=extents->index;
 853   /*
 854     Find the right side (minima).
 855   */
 856   for ( ; extents->index <= 255; extents->index++)
 857     if (extrema[extents->index] < 0)
 858       break;
 859   extents->right=extents->index-1;
 860   return(MagickTrue);
 861 }
 862 \f
 863 /*
 864 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 865 %                                                                             %
 866 %                                                                             %
 867 %                                                                             %
 868 +   D e r i v a t i v e H i s t o g r a m                                     %
 869 %                                                                             %
 870 %                                                                             %
 871 %                                                                             %
 872 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 873 %
 874 %  DerivativeHistogram() determines the derivative of the histogram using
 875 %  central differencing.
 876 %
 877 %  The format of the DerivativeHistogram method is:
 878 %
 879 %      DerivativeHistogram(const MagickRealType *histogram,
 880 %        MagickRealType *derivative)
 881 %
 882 %  A description of each parameter follows.
 883 %
 884 %    o histogram: Specifies an array of MagickRealTypes representing the number
 885 %      of pixels for each intensity of a particular color component.
 886 %
 887 %    o derivative: This array of MagickRealTypes is initialized by
 888 %      DerivativeHistogram to the derivative of the histogram using central
 889 %      differencing.
 890 %
 891 */
 892 static void DerivativeHistogram(const MagickRealType *histogram,
 893   MagickRealType *derivative)
 894 {
 895   register ssize_t
 896     i,
 897     n;
 898
 899   /*
 900     Compute endpoints using second order polynomial interpolation.
 901   */
 902   n=255;
 903   derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);
 904   derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);
 905   /*
 906     Compute derivative using central differencing.
 907   */
 908   for (i=1; i < n; i++)
 909     derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
 910   return;
 911 }
 912 \f
 913 /*
 914 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 915 %                                                                             %
 916 %                                                                             %
 917 %                                                                             %
 918 +  G e t I m a g e D y n a m i c T h r e s h o l d                            %
 919 %                                                                             %
 920 %                                                                             %
 921 %                                                                             %
 922 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 923 %
 924 %  GetImageDynamicThreshold() returns the dynamic threshold for an image.
 925 %
 926 %  The format of the GetImageDynamicThreshold method is:
 927 %
 928 %      MagickBooleanType GetImageDynamicThreshold(const Image *image,
 929 %        const double cluster_threshold,const double smooth_threshold,
 930 %        PixelInfo *pixel,ExceptionInfo *exception)
 931 %
 932 %  A description of each parameter follows.
 933 %
 934 %    o image: the image.
 935 %
 936 %    o cluster_threshold:  This MagickRealType represents the minimum number of
 937 %      pixels contained in a hexahedra before it can be considered valid
 938 %      (expressed as a percentage).
 939 %
 940 %    o smooth_threshold: the smoothing threshold eliminates noise in the second
 941 %      derivative of the histogram.  As the value is increased, you can expect a
 942 %      smoother second derivative.
 943 %
 944 %    o pixel: return the dynamic threshold here.
 945 %
 946 %    o exception: return any errors or warnings in this structure.
 947 %
 948 */
 949 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,
 950   const double cluster_threshold,const double smooth_threshold,
 951   PixelInfo *pixel,ExceptionInfo *exception)
 952 {
 953   Cluster
 954     *background,
 955     *cluster,
 956     *object,
 957     *head,
 958     *last_cluster,
 959     *next_cluster;
 960
 961   ExtentPacket
 962     blue,
 963     green,
 964     red;
 965
 966   MagickBooleanType
 967     proceed;
 968
 969   MagickRealType
 970     threshold;
 971
 972   register const Quantum
 973     *p;
 974
 975   register ssize_t
 976     i,
 977     x;
 978
 979   short
 980     *extrema[MaxDimension];
 981
 982   ssize_t
 983     count,
 984     *histogram[MaxDimension],
 985     y;
 986
 987   /*
 988     Allocate histogram and extrema.
 989   */
 990   assert(image != (Image *) NULL);
 991   assert(image->signature == MagickSignature);
 992   if (image->debug != MagickFalse)
 993     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 994   GetPixelInfo(image,pixel);
 995   for (i=0; i < MaxDimension; i++)
 996   {
 997     histogram[i]=(ssize_t *) AcquireQuantumMemory(256UL,sizeof(**histogram));
 998     extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
 999     if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1000       {
1001         for (i-- ; i >= 0; i--)
1002         {
1003           extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1004           histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1005         }
1006         (void) ThrowMagickException(exception,GetMagickModule(),
1007           ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1008         return(MagickFalse);
1009       }
1010   }
1011   /*
1012     Initialize histogram.
1013   */
1014   InitializeHistogram(image,histogram,exception);
1015   (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
1016     (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
1017   (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
1018     (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
1019   (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
1020     (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
1021   /*
1022     Form clusters.
1023   */
1024   cluster=(Cluster *) NULL;
1025   head=(Cluster *) NULL;
1026   (void) ResetMagickMemory(&red,0,sizeof(red));
1027   (void) ResetMagickMemory(&green,0,sizeof(green));
1028   (void) ResetMagickMemory(&blue,0,sizeof(blue));
1029   while (DefineRegion(extrema[Red],&red) != 0)
1030   {
1031     green.index=0;
1032     while (DefineRegion(extrema[Green],&green) != 0)
1033     {
1034       blue.index=0;
1035       while (DefineRegion(extrema[Blue],&blue) != 0)
1036       {
1037         /*
1038           Allocate a new class.
1039         */
1040         if (head != (Cluster *) NULL)
1041           {
1042             cluster->next=(Cluster *) AcquireMagickMemory(
1043               sizeof(*cluster->next));
1044             cluster=cluster->next;
1045           }
1046         else
1047           {
1048             cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1049             head=cluster;
1050           }
1051         if (cluster == (Cluster *) NULL)
1052           {
1053             (void) ThrowMagickException(exception,GetMagickModule(),
1054               ResourceLimitError,"MemoryAllocationFailed","`%s'",
1055               image->filename);
1056             return(MagickFalse);
1057           }
1058         /*
1059           Initialize a new class.
1060         */
1061         cluster->count=0;
1062         cluster->red=red;
1063         cluster->green=green;
1064         cluster->blue=blue;
1065         cluster->next=(Cluster *) NULL;
1066       }
1067     }
1068   }
1069   if (head == (Cluster *) NULL)
1070     {
1071       /*
1072         No classes were identified-- create one.
1073       */
1074       cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1075       if (cluster == (Cluster *) NULL)
1076         {
1077           (void) ThrowMagickException(exception,GetMagickModule(),
1078             ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1079           return(MagickFalse);
1080         }
1081       /*
1082         Initialize a new class.
1083       */
1084       cluster->count=0;
1085       cluster->red=red;
1086       cluster->green=green;
1087       cluster->blue=blue;
1088       cluster->next=(Cluster *) NULL;
1089       head=cluster;
1090     }
1091   /*
1092     Count the pixels for each cluster.
1093   */
1094   count=0;
1095   for (y=0; y < (ssize_t) image->rows; y++)
1096   {
1097     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1098     if (p == (const Quantum *) NULL)
1099       break;
1100     for (x=0; x < (ssize_t) image->columns; x++)
1101     {
1102       for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1103         if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >=
1104              (cluster->red.left-SafeMargin)) &&
1105             ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <=
1106              (cluster->red.right+SafeMargin)) &&
1107             ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >=
1108              (cluster->green.left-SafeMargin)) &&
1109             ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <=
1110              (cluster->green.right+SafeMargin)) &&
1111             ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >=
1112              (cluster->blue.left-SafeMargin)) &&
1113             ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <=
1114              (cluster->blue.right+SafeMargin)))
1115           {
1116             /*
1117               Count this pixel.
1118             */
1119             count++;
1120             cluster->red.center+=(MagickRealType) ScaleQuantumToChar(
1121               GetPixelRed(image,p));
1122             cluster->green.center+=(MagickRealType) ScaleQuantumToChar(
1123               GetPixelGreen(image,p));
1124             cluster->blue.center+=(MagickRealType) ScaleQuantumToChar(
1125               GetPixelBlue(image,p));
1126             cluster->count++;
1127             break;
1128           }
1129       p+=GetPixelChannels(image);
1130     }
1131     proceed=SetImageProgress(image,SegmentImageTag,(MagickOffsetType) y,
1132       2*image->rows);
1133     if (proceed == MagickFalse)
1134       break;
1135   }
1136   /*
1137     Remove clusters that do not meet minimum cluster threshold.
1138   */
1139   count=0;
1140   last_cluster=head;
1141   next_cluster=head;
1142   for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1143   {
1144     next_cluster=cluster->next;
1145     if ((cluster->count > 0) &&
1146         (cluster->count >= (count*cluster_threshold/100.0)))
1147       {
1148         /*
1149           Initialize cluster.
1150         */
1151         cluster->id=count;
1152         cluster->red.center/=cluster->count;
1153         cluster->green.center/=cluster->count;
1154         cluster->blue.center/=cluster->count;
1155         count++;
1156         last_cluster=cluster;
1157         continue;
1158       }
1159     /*
1160       Delete cluster.
1161     */
1162     if (cluster == head)
1163       head=next_cluster;
1164     else
1165       last_cluster->next=next_cluster;
1166     cluster=(Cluster *) RelinquishMagickMemory(cluster);
1167   }
1168   object=head;
1169   background=head;
1170   if (count > 1)
1171     {
1172       object=head->next;
1173       for (cluster=object; cluster->next != (Cluster *) NULL; )
1174       {
1175         if (cluster->count < object->count)
1176           object=cluster;
1177         cluster=cluster->next;
1178       }
1179       background=head->next;
1180       for (cluster=background; cluster->next != (Cluster *) NULL; )
1181       {
1182         if (cluster->count > background->count)
1183           background=cluster;
1184         cluster=cluster->next;
1185       }
1186     }
1187   threshold=(background->red.center+object->red.center)/2.0;
1188   pixel->red=(MagickRealType) ScaleCharToQuantum((unsigned char)
1189     (threshold+0.5));
1190   threshold=(background->green.center+object->green.center)/2.0;
1191   pixel->green=(MagickRealType) ScaleCharToQuantum((unsigned char)
1192     (threshold+0.5));
1193   threshold=(background->blue.center+object->blue.center)/2.0;
1194   pixel->blue=(MagickRealType) ScaleCharToQuantum((unsigned char)
1195     (threshold+0.5));
1196   /*
1197     Relinquish resources.
1198   */
1199   for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1200   {
1201     next_cluster=cluster->next;
1202     cluster=(Cluster *) RelinquishMagickMemory(cluster);
1203   }
1204   for (i=0; i < MaxDimension; i++)
1205   {
1206     extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1207     histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1208   }
1209   return(MagickTrue);
1210 }
1211 \f
1212 /*
1213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1214 %                                                                             %
1215 %                                                                             %
1216 %                                                                             %
1217 +  I n i t i a l i z e H i s t o g r a m                                      %
1218 %                                                                             %
1219 %                                                                             %
1220 %                                                                             %
1221 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1222 %
1223 %  InitializeHistogram() computes the histogram for an image.
1224 %
1225 %  The format of the InitializeHistogram method is:
1226 %
1227 %      InitializeHistogram(const Image *image,ssize_t **histogram)
1228 %
1229 %  A description of each parameter follows.
1230 %
1231 %    o image: Specifies a pointer to an Image structure;  returned from
1232 %      ReadImage.
1233 %
1234 %    o histogram: Specifies an array of integers representing the number
1235 %      of pixels for each intensity of a particular color component.
1236 %
1237 */
1238 static void InitializeHistogram(const Image *image,ssize_t **histogram,
1239   ExceptionInfo *exception)
1240 {
1241   register const Quantum
1242     *p;
1243
1244   register ssize_t
1245     i,
1246     x;
1247
1248   ssize_t
1249     y;
1250
1251   /*
1252     Initialize histogram.
1253   */
1254   for (i=0; i <= 255; i++)
1255   {
1256     histogram[Red][i]=0;
1257     histogram[Green][i]=0;
1258     histogram[Blue][i]=0;
1259   }
1260   for (y=0; y < (ssize_t) image->rows; y++)
1261   {
1262     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1263     if (p == (const Quantum *) NULL)
1264       break;
1265     for (x=0; x < (ssize_t) image->columns; x++)
1266     {
1267       histogram[Red][(ssize_t) ScaleQuantumToChar(GetPixelRed(image,p))]++;
1268       histogram[Green][(ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p))]++;
1269       histogram[Blue][(ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p))]++;
1270       p+=GetPixelChannels(image);
1271     }
1272   }
1273 }
1274 \f
1275 /*
1276 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1277 %                                                                             %
1278 %                                                                             %
1279 %                                                                             %
1280 +   I n i t i a l i z e I n t e r v a l T r e e                               %
1281 %                                                                             %
1282 %                                                                             %
1283 %                                                                             %
1284 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1285 %
1286 %  InitializeIntervalTree() initializes an interval tree from the lists of
1287 %  zero crossings.
1288 %
1289 %  The format of the InitializeIntervalTree method is:
1290 %
1291 %      InitializeIntervalTree(IntervalTree **list,ssize_t *number_nodes,
1292 %        IntervalTree *node)
1293 %
1294 %  A description of each parameter follows.
1295 %
1296 %    o zero_crossing: Specifies an array of structures of type ZeroCrossing.
1297 %
1298 %    o number_crossings: This size_t specifies the number of elements
1299 %      in the zero_crossing array.
1300 %
1301 */
1302
1303 static void InitializeList(IntervalTree **list,ssize_t *number_nodes,
1304   IntervalTree *node)
1305 {
1306   if (node == (IntervalTree *) NULL)
1307     return;
1308   if (node->child == (IntervalTree *) NULL)
1309     list[(*number_nodes)++]=node;
1310   InitializeList(list,number_nodes,node->sibling);
1311   InitializeList(list,number_nodes,node->child);
1312 }
1313
1314 static void MeanStability(IntervalTree *node)
1315 {
1316   register IntervalTree
1317     *child;
1318
1319   if (node == (IntervalTree *) NULL)
1320     return;
1321   node->mean_stability=0.0;
1322   child=node->child;
1323   if (child != (IntervalTree *) NULL)
1324     {
1325       register ssize_t
1326         count;
1327
1328       register MagickRealType
1329         sum;
1330
1331       sum=0.0;
1332       count=0;
1333       for ( ; child != (IntervalTree *) NULL; child=child->sibling)
1334       {
1335         sum+=child->stability;
1336         count++;
1337       }
1338       node->mean_stability=sum/(MagickRealType) count;
1339     }
1340   MeanStability(node->sibling);
1341   MeanStability(node->child);
1342 }
1343
1344 static void Stability(IntervalTree *node)
1345 {
1346   if (node == (IntervalTree *) NULL)
1347     return;
1348   if (node->child == (IntervalTree *) NULL)
1349     node->stability=0.0;
1350   else
1351     node->stability=node->tau-(node->child)->tau;
1352   Stability(node->sibling);
1353   Stability(node->child);
1354 }
1355
1356 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing,
1357   const size_t number_crossings)
1358 {
1359   IntervalTree
1360     *head,
1361     **list,
1362     *node,
1363     *root;
1364
1365   register ssize_t
1366     i;
1367
1368   ssize_t
1369     j,
1370     k,
1371     left,
1372     number_nodes;
1373
1374   /*
1375     Allocate interval tree.
1376   */
1377   list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1378     sizeof(*list));
1379   if (list == (IntervalTree **) NULL)
1380     return((IntervalTree *) NULL);
1381   /*
1382     The root is the entire histogram.
1383   */
1384   root=(IntervalTree *) AcquireMagickMemory(sizeof(*root));
1385   root->child=(IntervalTree *) NULL;
1386   root->sibling=(IntervalTree *) NULL;
1387   root->tau=0.0;
1388   root->left=0;
1389   root->right=255;
1390   for (i=(-1); i < (ssize_t) number_crossings; i++)
1391   {
1392     /*
1393       Initialize list with all nodes with no children.
1394     */
1395     number_nodes=0;
1396     InitializeList(list,&number_nodes,root);
1397     /*
1398       Split list.
1399     */
1400     for (j=0; j < number_nodes; j++)
1401     {
1402       head=list[j];
1403       left=head->left;
1404       node=head;
1405       for (k=head->left+1; k < head->right; k++)
1406       {
1407         if (zero_crossing[i+1].crossings[k] != 0)
1408           {
1409             if (node == head)
1410               {
1411                 node->child=(IntervalTree *) AcquireMagickMemory(
1412                   sizeof(*node->child));
1413                 node=node->child;
1414               }
1415             else
1416               {
1417                 node->sibling=(IntervalTree *) AcquireMagickMemory(
1418                   sizeof(*node->sibling));
1419                 node=node->sibling;
1420               }
1421             node->tau=zero_crossing[i+1].tau;
1422             node->child=(IntervalTree *) NULL;
1423             node->sibling=(IntervalTree *) NULL;
1424             node->left=left;
1425             node->right=k;
1426             left=k;
1427           }
1428         }
1429       if (left != head->left)
1430         {
1431           node->sibling=(IntervalTree *) AcquireMagickMemory(
1432             sizeof(*node->sibling));
1433           node=node->sibling;
1434           node->tau=zero_crossing[i+1].tau;
1435           node->child=(IntervalTree *) NULL;
1436           node->sibling=(IntervalTree *) NULL;
1437           node->left=left;
1438           node->right=head->right;
1439         }
1440     }
1441   }
1442   /*
1443     Determine the stability: difference between a nodes tau and its child.
1444   */
1445   Stability(root->child);
1446   MeanStability(root->child);
1447   list=(IntervalTree **) RelinquishMagickMemory(list);
1448   return(root);
1449 }
1450 \f
1451 /*
1452 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1453 %                                                                             %
1454 %                                                                             %
1455 %                                                                             %
1456 +   O p t i m a l T a u                                                       %
1457 %                                                                             %
1458 %                                                                             %
1459 %                                                                             %
1460 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1461 %
1462 %  OptimalTau() finds the optimal tau for each band of the histogram.
1463 %
1464 %  The format of the OptimalTau method is:
1465 %
1466 %    MagickRealType OptimalTau(const ssize_t *histogram,const double max_tau,
1467 %      const double min_tau,const double delta_tau,
1468 %      const double smooth_threshold,short *extrema)
1469 %
1470 %  A description of each parameter follows.
1471 %
1472 %    o histogram: Specifies an array of integers representing the number
1473 %      of pixels for each intensity of a particular color component.
1474 %
1475 %    o extrema:  Specifies a pointer to an array of integers.  They
1476 %      represent the peaks and valleys of the histogram for each color
1477 %      component.
1478 %
1479 */
1480
1481 static void ActiveNodes(IntervalTree **list,ssize_t *number_nodes,
1482   IntervalTree *node)
1483 {
1484   if (node == (IntervalTree *) NULL)
1485     return;
1486   if (node->stability >= node->mean_stability)
1487     {
1488       list[(*number_nodes)++]=node;
1489       ActiveNodes(list,number_nodes,node->sibling);
1490     }
1491   else
1492     {
1493       ActiveNodes(list,number_nodes,node->sibling);
1494       ActiveNodes(list,number_nodes,node->child);
1495     }
1496 }
1497
1498 static void FreeNodes(IntervalTree *node)
1499 {
1500   if (node == (IntervalTree *) NULL)
1501     return;
1502   FreeNodes(node->sibling);
1503   FreeNodes(node->child);
1504   node=(IntervalTree *) RelinquishMagickMemory(node);
1505 }
1506
1507 static MagickRealType OptimalTau(const ssize_t *histogram,const double max_tau,
1508   const double min_tau,const double delta_tau,const double smooth_threshold,
1509   short *extrema)
1510 {
1511   IntervalTree
1512     **list,
1513     *node,
1514     *root;
1515
1516   MagickBooleanType
1517     peak;
1518
1519   MagickRealType
1520     average_tau,
1521     *derivative,
1522     *second_derivative,
1523     tau,
1524     value;
1525
1526   register ssize_t
1527     i,
1528     x;
1529
1530   size_t
1531     count,
1532     number_crossings;
1533
1534   ssize_t
1535     index,
1536     j,
1537     k,
1538     number_nodes;
1539
1540   ZeroCrossing
1541     *zero_crossing;
1542
1543   /*
1544     Allocate interval tree.
1545   */
1546   list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1547     sizeof(*list));
1548   if (list == (IntervalTree **) NULL)
1549     return(0.0);
1550   /*
1551     Allocate zero crossing list.
1552   */
1553   count=(size_t) ((max_tau-min_tau)/delta_tau)+2;
1554   zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count,
1555     sizeof(*zero_crossing));
1556   if (zero_crossing == (ZeroCrossing *) NULL)
1557     return(0.0);
1558   for (i=0; i < (ssize_t) count; i++)
1559     zero_crossing[i].tau=(-1.0);
1560   /*
1561     Initialize zero crossing list.
1562   */
1563   derivative=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*derivative));
1564   second_derivative=(MagickRealType *) AcquireQuantumMemory(256,
1565     sizeof(*second_derivative));
1566   if ((derivative == (MagickRealType *) NULL) ||
1567       (second_derivative == (MagickRealType *) NULL))
1568     ThrowFatalException(ResourceLimitFatalError,
1569       "UnableToAllocateDerivatives");
1570   i=0;
1571   for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
1572   {
1573     zero_crossing[i].tau=tau;
1574     ScaleSpace(histogram,tau,zero_crossing[i].histogram);
1575     DerivativeHistogram(zero_crossing[i].histogram,derivative);
1576     DerivativeHistogram(derivative,second_derivative);
1577     ZeroCrossHistogram(second_derivative,smooth_threshold,
1578       zero_crossing[i].crossings);
1579     i++;
1580   }
1581   /*
1582     Add an entry for the original histogram.
1583   */
1584   zero_crossing[i].tau=0.0;
1585   for (j=0; j <= 255; j++)
1586     zero_crossing[i].histogram[j]=(MagickRealType) histogram[j];
1587   DerivativeHistogram(zero_crossing[i].histogram,derivative);
1588   DerivativeHistogram(derivative,second_derivative);
1589   ZeroCrossHistogram(second_derivative,smooth_threshold,
1590     zero_crossing[i].crossings);
1591   number_crossings=(size_t) i;
1592   derivative=(MagickRealType *) RelinquishMagickMemory(derivative);
1593   second_derivative=(MagickRealType *)
1594     RelinquishMagickMemory(second_derivative);
1595   /*
1596     Ensure the scale-space fingerprints form lines in scale-space, not loops.
1597   */
1598   ConsolidateCrossings(zero_crossing,number_crossings);
1599   /*
1600     Force endpoints to be included in the interval.
1601   */
1602   for (i=0; i <= (ssize_t) number_crossings; i++)
1603   {
1604     for (j=0; j < 255; j++)
1605       if (zero_crossing[i].crossings[j] != 0)
1606         break;
1607     zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]);
1608     for (j=255; j > 0; j--)
1609       if (zero_crossing[i].crossings[j] != 0)
1610         break;
1611     zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1612   }
1613   /*
1614     Initialize interval tree.
1615   */
1616   root=InitializeIntervalTree(zero_crossing,number_crossings);
1617   if (root == (IntervalTree *) NULL)
1618     return(0.0);
1619   /*
1620     Find active nodes:  stability is greater (or equal) to the mean stability of
1621     its children.
1622   */
1623   number_nodes=0;
1624   ActiveNodes(list,&number_nodes,root->child);
1625   /*
1626     Initialize extrema.
1627   */
1628   for (i=0; i <= 255; i++)
1629     extrema[i]=0;
1630   for (i=0; i < number_nodes; i++)
1631   {
1632     /*
1633       Find this tau in zero crossings list.
1634     */
1635     k=0;
1636     node=list[i];
1637     for (j=0; j <= (ssize_t) number_crossings; j++)
1638       if (zero_crossing[j].tau == node->tau)
1639         k=j;
1640     /*
1641       Find the value of the peak.
1642     */
1643     peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1644       MagickFalse;
1645     index=node->left;
1646     value=zero_crossing[k].histogram[index];
1647     for (x=node->left; x <= node->right; x++)
1648     {
1649       if (peak != MagickFalse)
1650         {
1651           if (zero_crossing[k].histogram[x] > value)
1652             {
1653               value=zero_crossing[k].histogram[x];
1654               index=x;
1655             }
1656         }
1657       else
1658         if (zero_crossing[k].histogram[x] < value)
1659           {
1660             value=zero_crossing[k].histogram[x];
1661             index=x;
1662           }
1663     }
1664     for (x=node->left; x <= node->right; x++)
1665     {
1666       if (index == 0)
1667         index=256;
1668       if (peak != MagickFalse)
1669         extrema[x]=(short) index;
1670       else
1671         extrema[x]=(short) (-index);
1672     }
1673   }
1674   /*
1675     Determine the average tau.
1676   */
1677   average_tau=0.0;
1678   for (i=0; i < number_nodes; i++)
1679     average_tau+=list[i]->tau;
1680   average_tau/=(MagickRealType) number_nodes;
1681   /*
1682     Relinquish resources.
1683   */
1684   FreeNodes(root);
1685   zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1686   list=(IntervalTree **) RelinquishMagickMemory(list);
1687   return(average_tau);
1688 }
1689 \f
1690 /*
1691 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1692 %                                                                             %
1693 %                                                                             %
1694 %                                                                             %
1695 +   S c a l e S p a c e                                                       %
1696 %                                                                             %
1697 %                                                                             %
1698 %                                                                             %
1699 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1700 %
1701 %  ScaleSpace() performs a scale-space filter on the 1D histogram.
1702 %
1703 %  The format of the ScaleSpace method is:
1704 %
1705 %      ScaleSpace(const ssize_t *histogram,const MagickRealType tau,
1706 %        MagickRealType *scale_histogram)
1707 %
1708 %  A description of each parameter follows.
1709 %
1710 %    o histogram: Specifies an array of MagickRealTypes representing the number
1711 %      of pixels for each intensity of a particular color component.
1712 %
1713 */
1714
1715 static void ScaleSpace(const ssize_t *histogram,const MagickRealType tau,
1716   MagickRealType *scale_histogram)
1717 {
1718   MagickRealType
1719     alpha,
1720     beta,
1721     *gamma,
1722     sum;
1723
1724   register ssize_t
1725     u,
1726     x;
1727
1728   gamma=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*gamma));
1729   if (gamma == (MagickRealType *) NULL)
1730     ThrowFatalException(ResourceLimitFatalError,
1731       "UnableToAllocateGammaMap");
1732   alpha=1.0/(tau*sqrt(2.0*MagickPI));
1733   beta=(-1.0/(2.0*tau*tau));
1734   for (x=0; x <= 255; x++)
1735     gamma[x]=0.0;
1736   for (x=0; x <= 255; x++)
1737   {
1738     gamma[x]=exp((double) beta*x*x);
1739     if (gamma[x] < MagickEpsilon)
1740       break;
1741   }
1742   for (x=0; x <= 255; x++)
1743   {
1744     sum=0.0;
1745     for (u=0; u <= 255; u++)
1746       sum+=(MagickRealType) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1747     scale_histogram[x]=alpha*sum;
1748   }
1749   gamma=(MagickRealType *) RelinquishMagickMemory(gamma);
1750 }
1751 \f
1752 /*
1753 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1754 %                                                                             %
1755 %                                                                             %
1756 %                                                                             %
1757 %  S e g m e n t I m a g e                                                    %
1758 %                                                                             %
1759 %                                                                             %
1760 %                                                                             %
1761 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1762 %
1763 %  SegmentImage() segment an image by analyzing the histograms of the color
1764 %  components and identifying units that are homogeneous with the fuzzy
1765 %  C-means technique.
1766 %
1767 %  The format of the SegmentImage method is:
1768 %
1769 %      MagickBooleanType SegmentImage(Image *image,
1770 %        const ColorspaceType colorspace,const MagickBooleanType verbose,
1771 %        const double cluster_threshold,const double smooth_threshold)
1772 %
1773 %  A description of each parameter follows.
1774 %
1775 %    o image: the image.
1776 %
1777 %    o colorspace: Indicate the colorspace.
1778 %
1779 %    o verbose:  Set to MagickTrue to print detailed information about the
1780 %      identified classes.
1781 %
1782 %    o cluster_threshold:  This represents the minimum number of pixels
1783 %      contained in a hexahedra before it can be considered valid (expressed
1784 %      as a percentage).
1785 %
1786 %    o smooth_threshold: the smoothing threshold eliminates noise in the second
1787 %      derivative of the histogram.  As the value is increased, you can expect a
1788 %      smoother second derivative.
1789 %
1790 */
1791 MagickExport MagickBooleanType SegmentImage(Image *image,
1792   const ColorspaceType colorspace,const MagickBooleanType verbose,
1793   const double cluster_threshold,const double smooth_threshold)
1794 {
1795   MagickBooleanType
1796     status;
1797
1798   register ssize_t
1799     i;
1800
1801   short
1802     *extrema[MaxDimension];
1803
1804   ssize_t
1805     *histogram[MaxDimension];
1806
1807   /*
1808     Allocate histogram and extrema.
1809   */
1810   assert(image != (Image *) NULL);
1811   assert(image->signature == MagickSignature);
1812   if (image->debug != MagickFalse)
1813     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1814   for (i=0; i < MaxDimension; i++)
1815   {
1816     histogram[i]=(ssize_t *) AcquireQuantumMemory(256,sizeof(**histogram));
1817     extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1818     if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1819       {
1820         for (i-- ; i >= 0; i--)
1821         {
1822           extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1823           histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1824         }
1825         ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1826           image->filename)
1827       }
1828   }
1829   if (IsRGBColorspace(colorspace) == MagickFalse)
1830     (void) TransformImageColorspace(image,colorspace);
1831   /*
1832     Initialize histogram.
1833   */
1834   InitializeHistogram(image,histogram,&image->exception);
1835   (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,
1836     smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Red]);
1837   (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,
1838     smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Green]);
1839   (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,
1840     smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Blue]);
1841   /*
1842     Classify using the fuzzy c-Means technique.
1843   */
1844   status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose);
1845   if (IsRGBColorspace(colorspace) == MagickFalse)
1846     (void) TransformImageColorspace(image,colorspace);
1847   /*
1848     Relinquish resources.
1849   */
1850   for (i=0; i < MaxDimension; i++)
1851   {
1852     extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1853     histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1854   }
1855   return(status);
1856 }
1857 \f
1858 /*
1859 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1860 %                                                                             %
1861 %                                                                             %
1862 %                                                                             %
1863 +   Z e r o C r o s s H i s t o g r a m                                       %
1864 %                                                                             %
1865 %                                                                             %
1866 %                                                                             %
1867 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1868 %
1869 %  ZeroCrossHistogram() find the zero crossings in a histogram and marks
1870 %  directions as:  1 is negative to positive; 0 is zero crossing; and -1
1871 %  is positive to negative.
1872 %
1873 %  The format of the ZeroCrossHistogram method is:
1874 %
1875 %      ZeroCrossHistogram(MagickRealType *second_derivative,
1876 %        const MagickRealType smooth_threshold,short *crossings)
1877 %
1878 %  A description of each parameter follows.
1879 %
1880 %    o second_derivative: Specifies an array of MagickRealTypes representing the
1881 %      second derivative of the histogram of a particular color component.
1882 %
1883 %    o crossings:  This array of integers is initialized with
1884 %      -1, 0, or 1 representing the slope of the first derivative of the
1885 %      of a particular color component.
1886 %
1887 */
1888 static void ZeroCrossHistogram(MagickRealType *second_derivative,
1889   const MagickRealType smooth_threshold,short *crossings)
1890 {
1891   register ssize_t
1892     i;
1893
1894   ssize_t
1895     parity;
1896
1897   /*
1898     Merge low numbers to zero to help prevent noise.
1899   */
1900   for (i=0; i <= 255; i++)
1901     if ((second_derivative[i] < smooth_threshold) &&
1902         (second_derivative[i] >= -smooth_threshold))
1903       second_derivative[i]=0.0;
1904   /*
1905     Mark zero crossings.
1906   */
1907   parity=0;
1908   for (i=0; i <= 255; i++)
1909   {
1910     crossings[i]=0;
1911     if (second_derivative[i] < 0.0)
1912       {
1913         if (parity > 0)
1914           crossings[i]=(-1);
1915         parity=1;
1916       }
1917     else
1918       if (second_derivative[i] > 0.0)
1919         {
1920           if (parity < 0)
1921             crossings[i]=1;
1922           parity=(-1);
1923         }
1924   }
1925 }