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