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