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