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