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