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