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