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