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"
111#define MaxDimension 3
113#if defined(FastClassify)
114#define WeightingExponent 2.0
115#define SegmentPower(ratio) (ratio)
117#define WeightingExponent 2.5
118#define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0)));
193 OptimalTau(
const ssize_t *,
const double,
const double,
const double,
194 const double,
short *);
202 ScaleSpace(
const ssize_t *,
const double,
double *),
203 ZeroCrossHistogram(
double *,
const double,
short *);
246static MagickBooleanType Classify(
Image *image,
short **extrema,
247 const double cluster_threshold,
const double weighting_exponent,
250#define SegmentImageTag "Segment/Image"
251#define ThrowClassifyException(severity,tag,label) \
253 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster) \
255 next_cluster=cluster->next; \
256 cluster=(Cluster *) RelinquishMagickMemory(cluster); \
258 if (squares != (double *) NULL) \
261 free_squares=squares; \
262 free_squares=(double *) RelinquishMagickMemory(free_squares); \
264 ThrowBinaryException(severity,tag,label); \
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)
315 while (DefineRegion(extrema[Green],&green) != 0)
318 while (DefineRegion(extrema[Blue],&blue) != 0)
325 cluster->next=(
Cluster *) AcquireQuantumMemory(1,
326 sizeof(*cluster->next));
327 cluster=cluster->next;
331 cluster=(
Cluster *) AcquireQuantumMemory(1,
sizeof(*cluster));
334 if (cluster == (
Cluster *) NULL)
335 ThrowClassifyException(ResourceLimitError,
"MemoryAllocationFailed",
340 (void) memset(cluster,0,
sizeof(*cluster));
342 cluster->green=green;
352 cluster=(
Cluster *) AcquireQuantumMemory(1,
sizeof(*cluster));
353 if (cluster == (
Cluster *) NULL)
354 ThrowClassifyException(ResourceLimitError,
"MemoryAllocationFailed",
359 (void) memset(cluster,0,
sizeof(*cluster));
361 cluster->green=green;
371 image_view=AcquireVirtualCacheView(image,exception);
372 for (y=0; y < (ssize_t) image->rows; y++)
380 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
381 if (p == (
const Quantum *) NULL)
383 for (x=0; x < (ssize_t) image->columns; x++)
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)))
403 cluster->red.center+=pixel.red;
404 cluster->green.center+=pixel.green;
405 cluster->blue.center+=pixel.blue;
409 p+=(ptrdiff_t) GetPixelChannels(image);
411 if (image->progress_monitor != (MagickProgressMonitor) NULL)
416#if defined(MAGICKCORE_OPENMP_SUPPORT)
420 proceed=SetImageProgress(image,SegmentImageTag,progress,2*image->rows);
421 if (proceed == MagickFalse)
425 image_view=DestroyCacheView(image_view);
432 for (cluster=head; cluster != (
Cluster *) NULL; cluster=next_cluster)
434 next_cluster=cluster->next;
435 if ((cluster->count > 0) &&
436 (cluster->count >= (count*cluster_threshold/100.0)))
442 cluster->red.center/=cluster->count;
443 cluster->green.center/=cluster->count;
444 cluster->blue.center/=cluster->count;
446 last_cluster=cluster;
455 last_cluster->next=next_cluster;
456 cluster=(
Cluster *) RelinquishMagickMemory(cluster);
458 number_clusters=(size_t) count;
459 if (verbose != MagickFalse)
464 (void) FormatLocaleFile(stdout,
"Fuzzy C-means Statistics\n");
465 (void) FormatLocaleFile(stdout,
"===================\n\n");
466 (void) FormatLocaleFile(stdout,
"\tCluster Threshold = %g\n",(
double)
468 (void) FormatLocaleFile(stdout,
"\tWeighting Exponent = %g\n",(
double)
470 (void) FormatLocaleFile(stdout,
"\tTotal Number of Clusters = %.20g\n\n",
471 (
double) number_clusters);
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);
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)
488 (void) FormatLocaleFile(stdout,
"\n\nCluster #%.20g\n\n",(
double)
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);
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)
504 (void) FormatLocaleFile(stdout,
"\n\nCluster #%.20g\n\n",(
double)
506 (void) FormatLocaleFile(stdout,
"%g %g %g\n",(
double)
507 cluster->red.center,(
double) cluster->green.center,(
double)
508 cluster->blue.center);
510 (void) FormatLocaleFile(stdout,
"\n");
512 if (number_clusters > 256)
513 ThrowClassifyException(ImageError,
"TooManyClusters",image->filename);
517 squares=(
double *) AcquireQuantumMemory(513UL,
sizeof(*squares));
518 if (squares == (
double *) NULL)
519 ThrowClassifyException(ResourceLimitError,
"MemoryAllocationFailed",
522 for (i=(-255); i <= 255; i++)
523 squares[i]=(
double) i*(double) i;
527 if (AcquireImageColormap(image,number_clusters,exception) == MagickFalse)
528 ThrowClassifyException(ResourceLimitError,
"MemoryAllocationFailed",
531 for (cluster=head; cluster != (
Cluster *) NULL; cluster=cluster->next)
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));
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)
549 for (y=0; y < (ssize_t) image->rows; y++)
563 if (status == MagickFalse)
565 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
566 if (q == (Quantum *) NULL)
571 for (x=0; x < (ssize_t) image->columns; x++)
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)
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)))
592 SetPixelIndex(image,(Quantum) c->id,q);
613 for (j=0; j < (ssize_t) image->colors; j++)
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++)
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);
632 if ((sum != 0.0) && ((1.0/sum) > local_minima))
637 local_minima=1.0/sum;
638 SetPixelIndex(image,(Quantum) j,q);
642 q+=(ptrdiff_t) GetPixelChannels(image);
644 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
646 if (image->progress_monitor != (MagickProgressMonitor) NULL)
651#if defined(MAGICKCORE_OPENMP_SUPPORT)
655 proceed=SetImageProgress(image,SegmentImageTag,progress,2*image->rows);
656 if (proceed == MagickFalse)
660 image_view=DestroyCacheView(image_view);
661 status&=(MagickStatusType) SyncImage(image,exception);
665 for (cluster=head; cluster != (
Cluster *) NULL; cluster=next_cluster)
667 next_cluster=cluster->next;
668 cluster=(
Cluster *) RelinquishMagickMemory(cluster);
671 free_squares=squares;
672 free_squares=(
double *) RelinquishMagickMemory(free_squares);
703static void ConsolidateCrossings(
ZeroCrossing *zero_crossing,
704 const size_t number_crossings)
722 for (i=(ssize_t) number_crossings-1; i >= 0; i--)
723 for (j=0; j <= 255; j++)
725 if (zero_crossing[i].crossings[j] == 0)
732 for (k=j-1; k > 0; k--)
733 if (zero_crossing[i+1].crossings[k] != 0)
737 for (k=j+1; k < 255; k++)
738 if (zero_crossing[i+1].crossings[k] != 0)
740 right=MagickMin(k,255);
744 for (k=j-1; k > 0; k--)
745 if (zero_crossing[i].crossings[k] != 0)
753 if (zero_crossing[i+1].crossings[j] != 0)
756 for (l=k+1; l < center; l++)
757 if (zero_crossing[i+1].crossings[l] != 0)
759 if (((count % 2) == 0) && (center != k))
768 for (l=k+1; l < left; l++)
769 if (zero_crossing[i+1].crossings[l] != 0)
771 if (((count % 2) == 0) && (left != k))
780 for (l=k+1; l < right; l++)
781 if (zero_crossing[i+1].crossings[l] != 0)
783 if (((count % 2) == 0) && (right != k))
786 l=(ssize_t) zero_crossing[i].crossings[j];
787 zero_crossing[i].crossings[j]=0;
789 zero_crossing[i].crossings[correct]=(short) l;
820static ssize_t DefineRegion(
const short *extrema,
ExtentPacket *extents)
831 for ( ; extents->index <= 255; extents->index++)
832 if (extrema[extents->index] > 0)
834 if (extents->index > 255)
836 extents->left=extents->index;
840 for ( ; extents->index <= 255; extents->index++)
841 if (extrema[extents->index] < 0)
843 extents->right=extents->index-1;
876static void DerivativeHistogram(
const double *histogram,
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]);
892 for (i=1; i < n; i++)
893 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
933MagickExport MagickBooleanType GetImageDynamicThreshold(
const Image *image,
934 const double cluster_threshold,
const double smooth_threshold,
964 *extrema[MaxDimension];
968 *histogram[MaxDimension],
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++)
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))
985 for (i-- ; i >= 0; i--)
987 extrema[i]=(
short *) RelinquishMagickMemory(extrema[i]);
988 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
990 (void) ThrowMagickException(exception,GetMagickModule(),
991 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
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]);
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)
1016 while (DefineRegion(extrema[Green],&green) != 0)
1019 while (DefineRegion(extrema[Blue],&blue) != 0)
1026 cluster->next=(
Cluster *) AcquireQuantumMemory(1,
1027 sizeof(*cluster->next));
1028 cluster=cluster->next;
1032 cluster=(
Cluster *) AcquireQuantumMemory(1,
sizeof(*cluster));
1035 if (cluster == (
Cluster *) NULL)
1037 (void) ThrowMagickException(exception,GetMagickModule(),
1038 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1040 return(MagickFalse);
1047 cluster->green=green;
1049 cluster->next=(
Cluster *) NULL;
1058 cluster=(
Cluster *) AcquireQuantumMemory(1,
sizeof(*cluster));
1059 if (cluster == (
Cluster *) NULL)
1061 (void) ThrowMagickException(exception,GetMagickModule(),
1062 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1063 return(MagickFalse);
1070 cluster->green=green;
1072 cluster->next=(
Cluster *) NULL;
1079 for (y=0; y < (ssize_t) image->rows; y++)
1081 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1082 if (p == (
const Quantum *) NULL)
1084 for (x=0; x < (ssize_t) image->columns; x++)
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)))
1106 cluster->red.center+=r;
1107 cluster->green.center+=g;
1108 cluster->blue.center+=b;
1112 p+=(ptrdiff_t) GetPixelChannels(image);
1114 proceed=SetImageProgress(image,SegmentImageTag,(MagickOffsetType) y,
1116 if (proceed == MagickFalse)
1125 for (cluster=head; cluster != (
Cluster *) NULL; cluster=next_cluster)
1127 next_cluster=cluster->next;
1128 if ((cluster->count > 0) &&
1129 (cluster->count >= (count*cluster_threshold/100.0)))
1135 cluster->red.center/=cluster->count;
1136 cluster->green.center/=cluster->count;
1137 cluster->blue.center/=cluster->count;
1139 last_cluster=cluster;
1145 if (cluster == head)
1148 last_cluster->next=next_cluster;
1149 cluster=(
Cluster *) RelinquishMagickMemory(cluster);
1153 if ((count > 1) && (head != (
Cluster *) NULL) &&
1154 (head->next != (
Cluster *) NULL))
1157 for (cluster=
object; cluster->next != (
Cluster *) NULL; )
1159 if (cluster->count < object->count)
1161 cluster=cluster->next;
1163 background=head->next;
1164 for (cluster=background; cluster->next != (
Cluster *) NULL; )
1166 if (cluster->count > background->count)
1168 cluster=cluster->next;
1171 if (background != (
Cluster *) NULL)
1173 threshold=(background->red.center+
object->red.center)/2.0;
1174 pixel->red=(double) ScaleCharToQuantum((
unsigned char)
1176 threshold=(background->green.center+
object->green.center)/2.0;
1177 pixel->green=(double) ScaleCharToQuantum((
unsigned char)
1179 threshold=(background->blue.center+
object->blue.center)/2.0;
1180 pixel->blue=(double) ScaleCharToQuantum((
unsigned char)
1186 for (cluster=head; cluster != (
Cluster *) NULL; cluster=next_cluster)
1188 next_cluster=cluster->next;
1189 cluster=(
Cluster *) RelinquishMagickMemory(cluster);
1191 for (i=0; i < MaxDimension; i++)
1193 extrema[i]=(
short *) RelinquishMagickMemory(extrema[i]);
1194 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1225static void InitializeHistogram(
const Image *image,ssize_t **histogram,
1241 for (i=0; i <= 255; i++)
1243 histogram[Red][i]=0;
1244 histogram[Green][i]=0;
1245 histogram[Blue][i]=0;
1247 for (y=0; y < (ssize_t) image->rows; y++)
1249 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1250 if (p == (
const Quantum *) NULL)
1252 for (x=0; x < (ssize_t) image->columns; x++)
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+=(ptrdiff_t) GetPixelChannels(image);
1290static void InitializeList(
IntervalTree **list,ssize_t *number_nodes,
1296 list[(*number_nodes)++]=node;
1297 InitializeList(list,number_nodes,node->sibling);
1298 InitializeList(list,number_nodes,node->child);
1308 node->mean_stability=0.0;
1320 for ( ; child != (
IntervalTree *) NULL; child=child->sibling)
1322 sum+=child->stability;
1325 node->mean_stability=sum/(double) count;
1327 MeanStability(node->sibling);
1328 MeanStability(node->child);
1336 node->stability=0.0;
1338 node->stability=node->tau-(node->child)->tau;
1339 Stability(node->sibling);
1340 Stability(node->child);
1344 const size_t number_crossings)
1364 list=(
IntervalTree **) AcquireQuantumMemory((
size_t) TreeLength,
1371 root=(
IntervalTree *) AcquireCriticalMemory(
sizeof(*root));
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++)
1386 InitializeList(list,&number_nodes,root);
1390 for (j=0; j < number_nodes; j++)
1395 for (k=head->left+1; k < head->right; k++)
1397 if (zero_crossing[i+1].crossings[k] != 0)
1402 sizeof(*node->child));
1408 sizeof(*node->sibling));
1417 node->tau=zero_crossing[i+1].tau;
1425 if (left != head->left)
1428 sizeof(*node->sibling));
1436 node->tau=zero_crossing[i+1].tau;
1440 node->right=head->right;
1447 Stability(root->child);
1448 MeanStability(root->child);
1483static void ActiveNodes(
IntervalTree **list,ssize_t *number_nodes,
1488 if (node->stability >= node->mean_stability)
1490 list[(*number_nodes)++]=node;
1491 ActiveNodes(list,number_nodes,node->sibling);
1495 ActiveNodes(list,number_nodes,node->sibling);
1496 ActiveNodes(list,number_nodes,node->child);
1504 FreeNodes(node->sibling);
1505 FreeNodes(node->child);
1509static double OptimalTau(
const ssize_t *histogram,
const double max_tau,
1510 const double min_tau,
const double delta_tau,
const double smooth_threshold,
1548 list=(
IntervalTree **) AcquireQuantumMemory((
size_t) TreeLength,
1555 count=(size_t) ((max_tau-min_tau)/delta_tau)+2;
1556 zero_crossing=(
ZeroCrossing *) AcquireQuantumMemory((
size_t) count,
1557 sizeof(*zero_crossing));
1563 for (i=0; i < (ssize_t) count; i++)
1564 zero_crossing[i].tau=(-1.0);
1568 derivative=(
double *) AcquireCriticalMemory(256*
sizeof(*derivative));
1569 second_derivative=(
double *) AcquireCriticalMemory(256*
1570 sizeof(*second_derivative));
1572 for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
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);
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);
1598 ConsolidateCrossings(zero_crossing,number_crossings);
1602 for (i=0; i <= (ssize_t) number_crossings; i++)
1604 for (j=0; j < 255; j++)
1605 if (zero_crossing[i].crossings[j] != 0)
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)
1611 zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1616 root=InitializeIntervalTree(zero_crossing,number_crossings);
1619 zero_crossing=(
ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1628 ActiveNodes(list,&number_nodes,root->child);
1632 for (i=0; i <= 255; i++)
1634 for (i=0; i < number_nodes; i++)
1641 for (j=0; j <= (ssize_t) number_crossings; j++)
1642 if (zero_crossing[j].tau == node->tau)
1647 peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1650 value=zero_crossing[k].histogram[index];
1651 for (x=node->left; x <= node->right; x++)
1653 if (peak != MagickFalse)
1655 if (zero_crossing[k].histogram[x] > value)
1657 value=zero_crossing[k].histogram[x];
1662 if (zero_crossing[k].histogram[x] < value)
1664 value=zero_crossing[k].histogram[x];
1668 for (x=node->left; x <= node->right; x++)
1672 if (peak != MagickFalse)
1673 extrema[x]=(short) index;
1675 extrema[x]=(short) (-index);
1682 for (i=0; i < number_nodes; i++)
1683 average_tau+=list[i]->tau;
1684 average_tau*=PerceptibleReciprocal((
double) number_nodes);
1689 zero_crossing=(
ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1691 return(average_tau);
1718static void ScaleSpace(
const ssize_t *histogram,
const double tau,
1719 double *scale_histogram)
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++)
1738 for (x=0; x <= 255; x++)
1740 gamma[x]=exp((
double) beta*x*x);
1741 if (gamma[x] < MagickEpsilon)
1744 for (x=0; x <= 255; x++)
1747 for (u=0; u <= 255; u++)
1748 sum+=(
double) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1749 scale_histogram[x]=alpha*sum;
1751 gamma=(
double *) RelinquishMagickMemory(gamma);
1796MagickExport MagickBooleanType SegmentImage(
Image *image,
1797 const ColorspaceType colorspace,
const MagickBooleanType verbose,
1798 const double cluster_threshold,
const double smooth_threshold,
1802 previous_colorspace;
1811 *extrema[MaxDimension];
1814 *histogram[MaxDimension];
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++)
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))
1829 for (i-- ; i >= 0; i--)
1831 extrema[i]=(
short *) RelinquishMagickMemory(extrema[i]);
1832 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1834 ThrowBinaryException(ResourceLimitError,
"MemoryAllocationFailed",
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]);
1853 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose,
1855 (void) TransformImageColorspace(image,previous_colorspace,exception);
1859 for (i=0; i < MaxDimension; i++)
1861 extrema[i]=(
short *) RelinquishMagickMemory(extrema[i]);
1862 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1897static void ZeroCrossHistogram(
double *second_derivative,
1898 const double smooth_threshold,
short *crossings)
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;
1917 for (i=0; i <= 255; i++)
1920 if (second_derivative[i] < 0.0)
1927 if (second_derivative[i] > 0.0)