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-2021 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,const double weighting_exponent,
225 % const MagickBooleanType verbose,ExceptionInfo *exception)
226 %
227 % A description of each parameter follows.
228 %
229 % o image: the image.
230 %
231 % o extrema: Specifies a pointer to an array of integers. They
232 % represent the peaks and valleys of the histogram for each color
233 % component.
234 %
235 % o cluster_threshold: This double represents the minimum number of
236 % pixels contained in a hexahedra before it can be considered valid
237 % (expressed as a percentage).
238 %
239 % o weighting_exponent: Specifies the membership weighting exponent.
240 %
241 % o verbose: A value greater than zero prints detailed information about
242 % the identified classes.
243 %
244 % o exception: return any errors or warnings in this structure.
245 %
246 */
Classify(Image * image,short ** extrema,const double cluster_threshold,const double weighting_exponent,const MagickBooleanType verbose,ExceptionInfo * exception)247 static MagickBooleanType Classify(Image *image,short **extrema,
248 const double cluster_threshold,const double weighting_exponent,
249 const MagickBooleanType verbose,ExceptionInfo *exception)
250 {
251 #define SegmentImageTag "Segment/Image"
252 #define ThrowClassifyException(severity,tag,label) \
253 {\
254 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster) \
255 { \
256 next_cluster=cluster->next; \
257 cluster=(Cluster *) RelinquishMagickMemory(cluster); \
258 } \
259 if (squares != (double *) NULL) \
260 { \
261 squares-=255; \
262 free_squares=squares; \
263 free_squares=(double *) RelinquishMagickMemory(free_squares); \
264 } \
265 ThrowBinaryException(severity,tag,label); \
266 }
267
268 CacheView
269 *image_view;
270
271 Cluster
272 *cluster,
273 *head,
274 *last_cluster,
275 *next_cluster;
276
277 double
278 *free_squares;
279
280 ExtentPacket
281 blue,
282 green,
283 red;
284
285 MagickOffsetType
286 progress;
287
288 MagickStatusType
289 status;
290
291 ssize_t
292 i;
293
294 double
295 *squares;
296
297 size_t
298 number_clusters;
299
300 ssize_t
301 count,
302 y;
303
304 /*
305 Form clusters.
306 */
307 cluster=(Cluster *) NULL;
308 head=(Cluster *) NULL;
309 squares=(double *) NULL;
310 (void) memset(&red,0,sizeof(red));
311 (void) memset(&green,0,sizeof(green));
312 (void) memset(&blue,0,sizeof(blue));
313 while (DefineRegion(extrema[Red],&red) != 0)
314 {
315 green.index=0;
316 while (DefineRegion(extrema[Green],&green) != 0)
317 {
318 blue.index=0;
319 while (DefineRegion(extrema[Blue],&blue) != 0)
320 {
321 /*
322 Allocate a new class.
323 */
324 if (head != (Cluster *) NULL)
325 {
326 cluster->next=(Cluster *) AcquireQuantumMemory(1,
327 sizeof(*cluster->next));
328 cluster=cluster->next;
329 }
330 else
331 {
332 cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));
333 head=cluster;
334 }
335 if (cluster == (Cluster *) NULL)
336 ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",
337 image->filename);
338 /*
339 Initialize a new class.
340 */
341 (void) memset(cluster,0,sizeof(*cluster));
342 cluster->red=red;
343 cluster->green=green;
344 cluster->blue=blue;
345 }
346 }
347 }
348 if (head == (Cluster *) NULL)
349 {
350 /*
351 No classes were identified-- create one.
352 */
353 cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));
354 if (cluster == (Cluster *) NULL)
355 ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",
356 image->filename);
357 /*
358 Initialize a new class.
359 */
360 (void) memset(cluster,0,sizeof(*cluster));
361 cluster->red=red;
362 cluster->green=green;
363 cluster->blue=blue;
364 head=cluster;
365 }
366 /*
367 Count the pixels for each cluster.
368 */
369 status=MagickTrue;
370 count=0;
371 progress=0;
372 image_view=AcquireVirtualCacheView(image,exception);
373 for (y=0; y < (ssize_t) image->rows; y++)
374 {
375 const Quantum
376 *p;
377
378 ssize_t
379 x;
380
381 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
382 if (p == (const Quantum *) NULL)
383 break;
384 for (x=0; x < (ssize_t) image->columns; x++)
385 {
386 PixelInfo
387 pixel;
388
389 pixel.red=(double) ScaleQuantumToChar(GetPixelRed(image,p));
390 pixel.green=(double) ScaleQuantumToChar(GetPixelGreen(image,p));
391 pixel.blue=(double) ScaleQuantumToChar(GetPixelBlue(image,p));
392 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
393 if ((pixel.red >= (double) (cluster->red.left-SafeMargin)) &&
394 (pixel.red <= (double) (cluster->red.right+SafeMargin)) &&
395 (pixel.green >= (double) (cluster->green.left-SafeMargin)) &&
396 (pixel.green <= (double) (cluster->green.right+SafeMargin)) &&
397 (pixel.blue >= (double) (cluster->blue.left-SafeMargin)) &&
398 (pixel.blue <= (double) (cluster->blue.right+SafeMargin)))
399 {
400 /*
401 Count this pixel.
402 */
403 count++;
404 cluster->red.center+=pixel.red;
405 cluster->green.center+=pixel.green;
406 cluster->blue.center+=pixel.blue;
407 cluster->count++;
408 break;
409 }
410 p+=GetPixelChannels(image);
411 }
412 if (image->progress_monitor != (MagickProgressMonitor) NULL)
413 {
414 MagickBooleanType
415 proceed;
416
417 #if defined(MAGICKCORE_OPENMP_SUPPORT)
418 #pragma omp atomic
419 #endif
420 progress++;
421 proceed=SetImageProgress(image,SegmentImageTag,progress,2*image->rows);
422 if (proceed == MagickFalse)
423 status=MagickFalse;
424 }
425 }
426 image_view=DestroyCacheView(image_view);
427 /*
428 Remove clusters that do not meet minimum cluster threshold.
429 */
430 count=0;
431 last_cluster=head;
432 next_cluster=head;
433 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
434 {
435 next_cluster=cluster->next;
436 if ((cluster->count > 0) &&
437 (cluster->count >= (count*cluster_threshold/100.0)))
438 {
439 /*
440 Initialize cluster.
441 */
442 cluster->id=count;
443 cluster->red.center/=cluster->count;
444 cluster->green.center/=cluster->count;
445 cluster->blue.center/=cluster->count;
446 count++;
447 last_cluster=cluster;
448 continue;
449 }
450 /*
451 Delete cluster.
452 */
453 if (cluster == head)
454 head=next_cluster;
455 else
456 last_cluster->next=next_cluster;
457 cluster=(Cluster *) RelinquishMagickMemory(cluster);
458 }
459 number_clusters=(size_t) count;
460 if (verbose != MagickFalse)
461 {
462 /*
463 Print cluster statistics.
464 */
465 (void) FormatLocaleFile(stdout,"Fuzzy C-means Statistics\n");
466 (void) FormatLocaleFile(stdout,"===================\n\n");
467 (void) FormatLocaleFile(stdout,"\tCluster Threshold = %g\n",(double)
468 cluster_threshold);
469 (void) FormatLocaleFile(stdout,"\tWeighting Exponent = %g\n",(double)
470 weighting_exponent);
471 (void) FormatLocaleFile(stdout,"\tTotal Number of Clusters = %.20g\n\n",
472 (double) number_clusters);
473 /*
474 Print the total number of points per cluster.
475 */
476 (void) FormatLocaleFile(stdout,"\n\nNumber of Vectors Per Cluster\n");
477 (void) FormatLocaleFile(stdout,"=============================\n\n");
478 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
479 (void) FormatLocaleFile(stdout,"Cluster #%.20g = %.20g\n",(double)
480 cluster->id,(double) cluster->count);
481 /*
482 Print the cluster extents.
483 */
484 (void) FormatLocaleFile(stdout,
485 "\n\n\nCluster Extents: (Vector Size: %d)\n",MaxDimension);
486 (void) FormatLocaleFile(stdout,"================");
487 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
488 {
489 (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
490 cluster->id);
491 (void) FormatLocaleFile(stdout,
492 "%.20g-%.20g %.20g-%.20g %.20g-%.20g\n",(double)
493 cluster->red.left,(double) cluster->red.right,(double)
494 cluster->green.left,(double) cluster->green.right,(double)
495 cluster->blue.left,(double) cluster->blue.right);
496 }
497 /*
498 Print the cluster center values.
499 */
500 (void) FormatLocaleFile(stdout,
501 "\n\n\nCluster Center Values: (Vector Size: %d)\n",MaxDimension);
502 (void) FormatLocaleFile(stdout,"=====================");
503 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
504 {
505 (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
506 cluster->id);
507 (void) FormatLocaleFile(stdout,"%g %g %g\n",(double)
508 cluster->red.center,(double) cluster->green.center,(double)
509 cluster->blue.center);
510 }
511 (void) FormatLocaleFile(stdout,"\n");
512 }
513 if (number_clusters > 256)
514 ThrowClassifyException(ImageError,"TooManyClusters",image->filename);
515 /*
516 Speed up distance calculations.
517 */
518 squares=(double *) AcquireQuantumMemory(513UL,sizeof(*squares));
519 if (squares == (double *) NULL)
520 ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",
521 image->filename);
522 squares+=255;
523 for (i=(-255); i <= 255; i++)
524 squares[i]=(double) i*(double) i;
525 /*
526 Allocate image colormap.
527 */
528 if (AcquireImageColormap(image,number_clusters,exception) == MagickFalse)
529 ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",
530 image->filename);
531 i=0;
532 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
533 {
534 image->colormap[i].red=(double) ScaleCharToQuantum((unsigned char)
535 (cluster->red.center+0.5));
536 image->colormap[i].green=(double) ScaleCharToQuantum((unsigned char)
537 (cluster->green.center+0.5));
538 image->colormap[i].blue=(double) ScaleCharToQuantum((unsigned char)
539 (cluster->blue.center+0.5));
540 i++;
541 }
542 /*
543 Do course grain classes.
544 */
545 image_view=AcquireAuthenticCacheView(image,exception);
546 #if defined(MAGICKCORE_OPENMP_SUPPORT)
547 #pragma omp parallel for schedule(static) shared(progress,status) \
548 magick_number_threads(image,image,image->rows,1)
549 #endif
550 for (y=0; y < (ssize_t) image->rows; y++)
551 {
552 Cluster
553 *cluster;
554
555 const PixelInfo
556 *magick_restrict p;
557
558 ssize_t
559 x;
560
561 Quantum
562 *magick_restrict q;
563
564 if (status == MagickFalse)
565 continue;
566 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
567 if (q == (Quantum *) NULL)
568 {
569 status=MagickFalse;
570 continue;
571 }
572 for (x=0; x < (ssize_t) image->columns; x++)
573 {
574 PixelInfo
575 pixel;
576
577 SetPixelIndex(image,(Quantum) 0,q);
578 pixel.red=(double) ScaleQuantumToChar(GetPixelRed(image,q));
579 pixel.green=(double) ScaleQuantumToChar(GetPixelGreen(image,q));
580 pixel.blue=(double) ScaleQuantumToChar(GetPixelBlue(image,q));
581 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
582 {
583 if ((pixel.red >= (double) (cluster->red.left-SafeMargin)) &&
584 (pixel.red <= (double) (cluster->red.right+SafeMargin)) &&
585 (pixel.green >= (double) (cluster->green.left-SafeMargin)) &&
586 (pixel.green <= (double) (cluster->green.right+SafeMargin)) &&
587 (pixel.blue >= (double) (cluster->blue.left-SafeMargin)) &&
588 (pixel.blue <= (double) (cluster->blue.right+SafeMargin)))
589 {
590 /*
591 Classify this pixel.
592 */
593 SetPixelIndex(image,(Quantum) cluster->id,q);
594 break;
595 }
596 }
597 if (cluster == (Cluster *) NULL)
598 {
599 double
600 distance_squared,
601 local_minima,
602 numerator,
603 ratio,
604 sum;
605
606 ssize_t
607 j,
608 k;
609
610 /*
611 Compute fuzzy membership.
612 */
613 local_minima=0.0;
614 for (j=0; j < (ssize_t) image->colors; j++)
615 {
616 sum=0.0;
617 p=image->colormap+j;
618 distance_squared=
619 squares[(ssize_t) (pixel.red-ScaleQuantumToChar(p->red))]+
620 squares[(ssize_t) (pixel.green-ScaleQuantumToChar(p->green))]+
621 squares[(ssize_t) (pixel.blue-ScaleQuantumToChar(p->blue))];
622 numerator=distance_squared;
623 for (k=0; k < (ssize_t) image->colors; k++)
624 {
625 p=image->colormap+k;
626 distance_squared=
627 squares[(ssize_t) (pixel.red-ScaleQuantumToChar(p->red))]+
628 squares[(ssize_t) (pixel.green-ScaleQuantumToChar(p->green))]+
629 squares[(ssize_t) (pixel.blue-ScaleQuantumToChar(p->blue))];
630 ratio=numerator/distance_squared;
631 sum+=SegmentPower(ratio);
632 }
633 if ((sum != 0.0) && ((1.0/sum) > local_minima))
634 {
635 /*
636 Classify this pixel.
637 */
638 local_minima=1.0/sum;
639 SetPixelIndex(image,(Quantum) j,q);
640 }
641 }
642 }
643 q+=GetPixelChannels(image);
644 }
645 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
646 status=MagickFalse;
647 if (image->progress_monitor != (MagickProgressMonitor) NULL)
648 {
649 MagickBooleanType
650 proceed;
651
652 #if defined(MAGICKCORE_OPENMP_SUPPORT)
653 #pragma omp atomic
654 #endif
655 progress++;
656 proceed=SetImageProgress(image,SegmentImageTag,progress,2*image->rows);
657 if (proceed == MagickFalse)
658 status=MagickFalse;
659 }
660 }
661 image_view=DestroyCacheView(image_view);
662 status&=SyncImage(image,exception);
663 /*
664 Relinquish resources.
665 */
666 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
667 {
668 next_cluster=cluster->next;
669 cluster=(Cluster *) RelinquishMagickMemory(cluster);
670 }
671 squares-=255;
672 free_squares=squares;
673 free_squares=(double *) RelinquishMagickMemory(free_squares);
674 return(MagickTrue);
675 }
676
677 /*
678 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
679 % %
680 % %
681 % %
682 + C o n s o l i d a t e C r o s s i n g s %
683 % %
684 % %
685 % %
686 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
687 %
688 % ConsolidateCrossings() guarantees that an even number of zero crossings
689 % always lie between two crossings.
690 %
691 % The format of the ConsolidateCrossings method is:
692 %
693 % ConsolidateCrossings(ZeroCrossing *zero_crossing,
694 % const size_t number_crossings)
695 %
696 % A description of each parameter follows.
697 %
698 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
699 %
700 % o number_crossings: This size_t specifies the number of elements
701 % in the zero_crossing array.
702 %
703 */
ConsolidateCrossings(ZeroCrossing * zero_crossing,const size_t number_crossings)704 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
705 const size_t number_crossings)
706 {
707 ssize_t
708 i,
709 j,
710 k,
711 l;
712
713 ssize_t
714 center,
715 correct,
716 count,
717 left,
718 right;
719
720 /*
721 Consolidate zero crossings.
722 */
723 for (i=(ssize_t) number_crossings-1; i >= 0; i--)
724 for (j=0; j <= 255; j++)
725 {
726 if (zero_crossing[i].crossings[j] == 0)
727 continue;
728 /*
729 Find the entry that is closest to j and still preserves the
730 property that there are an even number of crossings between
731 intervals.
732 */
733 for (k=j-1; k > 0; k--)
734 if (zero_crossing[i+1].crossings[k] != 0)
735 break;
736 left=MagickMax(k,0);
737 center=j;
738 for (k=j+1; k < 255; k++)
739 if (zero_crossing[i+1].crossings[k] != 0)
740 break;
741 right=MagickMin(k,255);
742 /*
743 K is the zero crossing just left of j.
744 */
745 for (k=j-1; k > 0; k--)
746 if (zero_crossing[i].crossings[k] != 0)
747 break;
748 if (k < 0)
749 k=0;
750 /*
751 Check center for an even number of crossings between k and j.
752 */
753 correct=(-1);
754 if (zero_crossing[i+1].crossings[j] != 0)
755 {
756 count=0;
757 for (l=k+1; l < center; l++)
758 if (zero_crossing[i+1].crossings[l] != 0)
759 count++;
760 if (((count % 2) == 0) && (center != k))
761 correct=center;
762 }
763 /*
764 Check left for an even number of crossings between k and j.
765 */
766 if (correct == -1)
767 {
768 count=0;
769 for (l=k+1; l < left; l++)
770 if (zero_crossing[i+1].crossings[l] != 0)
771 count++;
772 if (((count % 2) == 0) && (left != k))
773 correct=left;
774 }
775 /*
776 Check right for an even number of crossings between k and j.
777 */
778 if (correct == -1)
779 {
780 count=0;
781 for (l=k+1; l < right; l++)
782 if (zero_crossing[i+1].crossings[l] != 0)
783 count++;
784 if (((count % 2) == 0) && (right != k))
785 correct=right;
786 }
787 l=(ssize_t) zero_crossing[i].crossings[j];
788 zero_crossing[i].crossings[j]=0;
789 if (correct != -1)
790 zero_crossing[i].crossings[correct]=(short) l;
791 }
792 }
793
794 /*
795 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
796 % %
797 % %
798 % %
799 + D e f i n e R e g i o n %
800 % %
801 % %
802 % %
803 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
804 %
805 % DefineRegion() defines the left and right boundaries of a peak region.
806 %
807 % The format of the DefineRegion method is:
808 %
809 % ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
810 %
811 % A description of each parameter follows.
812 %
813 % o extrema: Specifies a pointer to an array of integers. They
814 % represent the peaks and valleys of the histogram for each color
815 % component.
816 %
817 % o extents: This pointer to an ExtentPacket represent the extends
818 % of a particular peak or valley of a color component.
819 %
820 */
DefineRegion(const short * extrema,ExtentPacket * extents)821 static ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
822 {
823 /*
824 Initialize to default values.
825 */
826 extents->left=0;
827 extents->center=0.0;
828 extents->right=255;
829 /*
830 Find the left side (maxima).
831 */
832 for ( ; extents->index <= 255; extents->index++)
833 if (extrema[extents->index] > 0)
834 break;
835 if (extents->index > 255)
836 return(MagickFalse); /* no left side - no region exists */
837 extents->left=extents->index;
838 /*
839 Find the right side (minima).
840 */
841 for ( ; extents->index <= 255; extents->index++)
842 if (extrema[extents->index] < 0)
843 break;
844 extents->right=extents->index-1;
845 return(MagickTrue);
846 }
847
848 /*
849 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
850 % %
851 % %
852 % %
853 + D e r i v a t i v e H i s t o g r a m %
854 % %
855 % %
856 % %
857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
858 %
859 % DerivativeHistogram() determines the derivative of the histogram using
860 % central differencing.
861 %
862 % The format of the DerivativeHistogram method is:
863 %
864 % DerivativeHistogram(const double *histogram,
865 % double *derivative)
866 %
867 % A description of each parameter follows.
868 %
869 % o histogram: Specifies an array of doubles representing the number
870 % of pixels for each intensity of a particular color component.
871 %
872 % o derivative: This array of doubles is initialized by
873 % DerivativeHistogram to the derivative of the histogram using central
874 % differencing.
875 %
876 */
DerivativeHistogram(const double * histogram,double * derivative)877 static void DerivativeHistogram(const double *histogram,
878 double *derivative)
879 {
880 ssize_t
881 i,
882 n;
883
884 /*
885 Compute endpoints using second order polynomial interpolation.
886 */
887 n=255;
888 derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);
889 derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);
890 /*
891 Compute derivative using central differencing.
892 */
893 for (i=1; i < n; i++)
894 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
895 return;
896 }
897
898 /*
899 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
900 % %
901 % %
902 % %
903 + G e t I m a g e D y n a m i c T h r e s h o l d %
904 % %
905 % %
906 % %
907 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
908 %
909 % GetImageDynamicThreshold() returns the dynamic threshold for an image.
910 %
911 % The format of the GetImageDynamicThreshold method is:
912 %
913 % MagickBooleanType GetImageDynamicThreshold(const Image *image,
914 % const double cluster_threshold,const double smooth_threshold,
915 % PixelInfo *pixel,ExceptionInfo *exception)
916 %
917 % A description of each parameter follows.
918 %
919 % o image: the image.
920 %
921 % o cluster_threshold: This double represents the minimum number of
922 % pixels contained in a hexahedra before it can be considered valid
923 % (expressed as a percentage).
924 %
925 % o smooth_threshold: the smoothing threshold eliminates noise in the second
926 % derivative of the histogram. As the value is increased, you can expect a
927 % smoother second derivative.
928 %
929 % o pixel: return the dynamic threshold here.
930 %
931 % o exception: return any errors or warnings in this structure.
932 %
933 */
GetImageDynamicThreshold(const Image * image,const double cluster_threshold,const double smooth_threshold,PixelInfo * pixel,ExceptionInfo * exception)934 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,
935 const double cluster_threshold,const double smooth_threshold,
936 PixelInfo *pixel,ExceptionInfo *exception)
937 {
938 Cluster
939 *background,
940 *cluster,
941 *object,
942 *head,
943 *last_cluster,
944 *next_cluster;
945
946 ExtentPacket
947 blue,
948 green,
949 red;
950
951 MagickBooleanType
952 proceed;
953
954 double
955 threshold;
956
957 const Quantum
958 *p;
959
960 ssize_t
961 i,
962 x;
963
964 short
965 *extrema[MaxDimension];
966
967 ssize_t
968 count,
969 *histogram[MaxDimension],
970 y;
971
972 /*
973 Allocate histogram and extrema.
974 */
975 assert(image != (Image *) NULL);
976 assert(image->signature == MagickCoreSignature);
977 if (image->debug != MagickFalse)
978 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
979 GetPixelInfo(image,pixel);
980 for (i=0; i < MaxDimension; i++)
981 {
982 histogram[i]=(ssize_t *) AcquireQuantumMemory(256UL,sizeof(**histogram));
983 extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
984 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
985 {
986 for (i-- ; i >= 0; i--)
987 {
988 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
989 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
990 }
991 (void) ThrowMagickException(exception,GetMagickModule(),
992 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
993 return(MagickFalse);
994 }
995 }
996 /*
997 Initialize histogram.
998 */
999 InitializeHistogram(image,histogram,exception);
1000 (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
1001 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
1002 (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
1003 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
1004 (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
1005 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
1006 /*
1007 Form clusters.
1008 */
1009 cluster=(Cluster *) NULL;
1010 head=(Cluster *) NULL;
1011 (void) memset(&red,0,sizeof(red));
1012 (void) memset(&green,0,sizeof(green));
1013 (void) memset(&blue,0,sizeof(blue));
1014 while (DefineRegion(extrema[Red],&red) != 0)
1015 {
1016 green.index=0;
1017 while (DefineRegion(extrema[Green],&green) != 0)
1018 {
1019 blue.index=0;
1020 while (DefineRegion(extrema[Blue],&blue) != 0)
1021 {
1022 /*
1023 Allocate a new class.
1024 */
1025 if (head != (Cluster *) NULL)
1026 {
1027 cluster->next=(Cluster *) AcquireQuantumMemory(1,
1028 sizeof(*cluster->next));
1029 cluster=cluster->next;
1030 }
1031 else
1032 {
1033 cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));
1034 head=cluster;
1035 }
1036 if (cluster == (Cluster *) NULL)
1037 {
1038 (void) ThrowMagickException(exception,GetMagickModule(),
1039 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1040 image->filename);
1041 return(MagickFalse);
1042 }
1043 /*
1044 Initialize a new class.
1045 */
1046 cluster->count=0;
1047 cluster->red=red;
1048 cluster->green=green;
1049 cluster->blue=blue;
1050 cluster->next=(Cluster *) NULL;
1051 }
1052 }
1053 }
1054 if (head == (Cluster *) NULL)
1055 {
1056 /*
1057 No classes were identified-- create one.
1058 */
1059 cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));
1060 if (cluster == (Cluster *) NULL)
1061 {
1062 (void) ThrowMagickException(exception,GetMagickModule(),
1063 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1064 return(MagickFalse);
1065 }
1066 /*
1067 Initialize a new class.
1068 */
1069 cluster->count=0;
1070 cluster->red=red;
1071 cluster->green=green;
1072 cluster->blue=blue;
1073 cluster->next=(Cluster *) NULL;
1074 head=cluster;
1075 }
1076 /*
1077 Count the pixels for each cluster.
1078 */
1079 count=0;
1080 for (y=0; y < (ssize_t) image->rows; y++)
1081 {
1082 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1083 if (p == (const Quantum *) NULL)
1084 break;
1085 for (x=0; x < (ssize_t) image->columns; x++)
1086 {
1087 PixelInfo
1088 pixel;
1089
1090 pixel.red=(double) ScaleQuantumToChar(GetPixelRed(image,p));
1091 pixel.green=(double) ScaleQuantumToChar(GetPixelGreen(image,p));
1092 pixel.blue=(double) ScaleQuantumToChar(GetPixelBlue(image,p));
1093 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1094 if ((pixel.red >= (double) (cluster->red.left-SafeMargin)) &&
1095 (pixel.red <= (double) (cluster->red.right+SafeMargin)) &&
1096 (pixel.green >= (double) (cluster->green.left-SafeMargin)) &&
1097 (pixel.green <= (double) (cluster->green.right+SafeMargin)) &&
1098 (pixel.blue >= (double) (cluster->blue.left-SafeMargin)) &&
1099 (pixel.blue <= (double) (cluster->blue.right+SafeMargin)))
1100 {
1101 /*
1102 Count this pixel.
1103 */
1104 count++;
1105 cluster->red.center+=pixel.red;
1106 cluster->green.center+=pixel.green;
1107 cluster->blue.center+=pixel.blue;
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 const Quantum
1227 *p;
1228
1229 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 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 ssize_t
1311 count;
1312
1313 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 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 *) AcquireQuantumMemory(1,
1400 sizeof(*node->child));
1401 node=node->child;
1402 }
1403 else
1404 {
1405 node->sibling=(IntervalTree *) AcquireQuantumMemory(1,
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 *) AcquireQuantumMemory(1,
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 double
1512 average_tau,
1513 *derivative,
1514 *second_derivative,
1515 tau,
1516 value;
1517
1518 IntervalTree
1519 **list,
1520 *node,
1521 *root;
1522
1523 MagickBooleanType
1524 peak;
1525
1526 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*=PerceptibleReciprocal((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 */
ScaleSpace(const ssize_t * histogram,const double tau,double * scale_histogram)1716 static void ScaleSpace(const ssize_t *histogram,const double tau,
1717 double *scale_histogram)
1718 {
1719 double
1720 alpha,
1721 beta,
1722 *gamma,
1723 sum;
1724
1725 ssize_t
1726 u,
1727 x;
1728
1729 gamma=(double *) AcquireQuantumMemory(256,sizeof(*gamma));
1730 if (gamma == (double *) NULL)
1731 ThrowFatalException(ResourceLimitFatalError,"UnableToAllocateGammaMap");
1732 alpha=PerceptibleReciprocal(tau*sqrt(2.0*MagickPI));
1733 beta=(-1.0*PerceptibleReciprocal(2.0*tau*tau));
1734 for (x=0; x <= 255; x++)
1735 gamma[x]=0.0;
1736 for (x=0; x <= 255; x++)
1737 {
1738 gamma[x]=exp((double) beta*x*x);
1739 if (gamma[x] < MagickEpsilon)
1740 break;
1741 }
1742 for (x=0; x <= 255; x++)
1743 {
1744 sum=0.0;
1745 for (u=0; u <= 255; u++)
1746 sum+=(double) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1747 scale_histogram[x]=alpha*sum;
1748 }
1749 gamma=(double *) RelinquishMagickMemory(gamma);
1750 }
1751
1752 /*
1753 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1754 % %
1755 % %
1756 % %
1757 % S e g m e n t I m a g e %
1758 % %
1759 % %
1760 % %
1761 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1762 %
1763 % SegmentImage() segment an image by analyzing the histograms of the color
1764 % components and identifying units that are homogeneous with the fuzzy
1765 % C-means technique.
1766 %
1767 % The format of the SegmentImage method is:
1768 %
1769 % MagickBooleanType SegmentImage(Image *image,
1770 % const ColorspaceType colorspace,const MagickBooleanType verbose,
1771 % const double cluster_threshold,const double smooth_threshold,
1772 % ExceptionInfo *exception)
1773 %
1774 % A description of each parameter follows.
1775 %
1776 % o image: the image.
1777 %
1778 % o colorspace: Indicate the colorspace.
1779 %
1780 % o verbose: Set to MagickTrue to print detailed information about the
1781 % identified classes.
1782 %
1783 % o cluster_threshold: This represents the minimum number of pixels
1784 % contained in a hexahedra before it can be considered valid (expressed
1785 % as a percentage).
1786 %
1787 % o smooth_threshold: the smoothing threshold eliminates noise in the second
1788 % derivative of the histogram. As the value is increased, you can expect a
1789 % smoother second derivative.
1790 %
1791 % o exception: return any errors or warnings in this structure.
1792 %
1793 */
SegmentImage(Image * image,const ColorspaceType colorspace,const MagickBooleanType verbose,const double cluster_threshold,const double smooth_threshold,ExceptionInfo * exception)1794 MagickExport MagickBooleanType SegmentImage(Image *image,
1795 const ColorspaceType colorspace,const MagickBooleanType verbose,
1796 const double cluster_threshold,const double smooth_threshold,
1797 ExceptionInfo *exception)
1798 {
1799 ColorspaceType
1800 previous_colorspace;
1801
1802 MagickBooleanType
1803 status;
1804
1805 ssize_t
1806 i;
1807
1808 short
1809 *extrema[MaxDimension];
1810
1811 ssize_t
1812 *histogram[MaxDimension];
1813
1814 /*
1815 Allocate histogram and extrema.
1816 */
1817 assert(image != (Image *) NULL);
1818 assert(image->signature == MagickCoreSignature);
1819 if (image->debug != MagickFalse)
1820 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1821 for (i=0; i < MaxDimension; i++)
1822 {
1823 histogram[i]=(ssize_t *) AcquireQuantumMemory(256,sizeof(**histogram));
1824 extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1825 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1826 {
1827 for (i-- ; i >= 0; i--)
1828 {
1829 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1830 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1831 }
1832 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1833 image->filename)
1834 }
1835 }
1836 /*
1837 Initialize histogram.
1838 */
1839 previous_colorspace=image->colorspace;
1840 (void) TransformImageColorspace(image,colorspace,exception);
1841 InitializeHistogram(image,histogram,exception);
1842 (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,smooth_threshold == 0.0 ?
1843 1.0 : smooth_threshold,extrema[Red]);
1844 (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,smooth_threshold == 0.0 ?
1845 1.0 : smooth_threshold,extrema[Green]);
1846 (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,smooth_threshold == 0.0 ?
1847 1.0 : smooth_threshold,extrema[Blue]);
1848 /*
1849 Classify using the fuzzy c-Means technique.
1850 */
1851 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose,
1852 exception);
1853 (void) TransformImageColorspace(image,previous_colorspace,exception);
1854 /*
1855 Relinquish resources.
1856 */
1857 for (i=0; i < MaxDimension; i++)
1858 {
1859 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1860 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1861 }
1862 return(status);
1863 }
1864
1865 /*
1866 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1867 % %
1868 % %
1869 % %
1870 + Z e r o C r o s s H i s t o g r a m %
1871 % %
1872 % %
1873 % %
1874 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1875 %
1876 % ZeroCrossHistogram() find the zero crossings in a histogram and marks
1877 % directions as: 1 is negative to positive; 0 is zero crossing; and -1
1878 % is positive to negative.
1879 %
1880 % The format of the ZeroCrossHistogram method is:
1881 %
1882 % ZeroCrossHistogram(double *second_derivative,
1883 % const double smooth_threshold,short *crossings)
1884 %
1885 % A description of each parameter follows.
1886 %
1887 % o second_derivative: Specifies an array of doubles representing the
1888 % second derivative of the histogram of a particular color component.
1889 %
1890 % o crossings: This array of integers is initialized with
1891 % -1, 0, or 1 representing the slope of the first derivative of the
1892 % of a particular color component.
1893 %
1894 */
ZeroCrossHistogram(double * second_derivative,const double smooth_threshold,short * crossings)1895 static void ZeroCrossHistogram(double *second_derivative,
1896 const double smooth_threshold,short *crossings)
1897 {
1898 ssize_t
1899 i;
1900
1901 ssize_t
1902 parity;
1903
1904 /*
1905 Merge low numbers to zero to help prevent noise.
1906 */
1907 for (i=0; i <= 255; i++)
1908 if ((second_derivative[i] < smooth_threshold) &&
1909 (second_derivative[i] >= -smooth_threshold))
1910 second_derivative[i]=0.0;
1911 /*
1912 Mark zero crossings.
1913 */
1914 parity=0;
1915 for (i=0; i <= 255; i++)
1916 {
1917 crossings[i]=0;
1918 if (second_derivative[i] < 0.0)
1919 {
1920 if (parity > 0)
1921 crossings[i]=(-1);
1922 parity=1;
1923 }
1924 else
1925 if (second_derivative[i] > 0.0)
1926 {
1927 if (parity < 0)
1928 crossings[i]=1;
1929 parity=(-1);
1930 }
1931 }
1932 }
1933