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