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