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