feature.c revision b6ecf44155ec6e79fac7d9c51ed5f56962f6eda4
1/* 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3% % 4% % 5% % 6% FFFFF EEEEE AAA TTTTT U U RRRR EEEEE % 7% F E A A T U U R R E % 8% FFF EEE AAAAA T U U RRRR EEE % 9% F E A A T U U R R E % 10% F EEEEE A A T UUU R R EEEEE % 11% % 12% % 13% MagickCore Image Feature Methods % 14% % 15% Software Design % 16% Cristy % 17% July 1992 % 18% % 19% % 20% Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization % 21% dedicated to making software imaging solutions freely available. % 22% % 23% You may not use this file except in compliance with the License. You may % 24% obtain a copy of the License at % 25% % 26% http://www.imagemagick.org/script/license.php % 27% % 28% Unless required by applicable law or agreed to in writing, software % 29% distributed under the License is distributed on an "AS IS" BASIS, % 30% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. % 31% See the License for the specific language governing permissions and % 32% limitations under the License. % 33% % 34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 35% 36% 37% 38*/ 39 40/* 41 Include declarations. 42*/ 43#include "MagickCore/studio.h" 44#include "MagickCore/property.h" 45#include "MagickCore/animate.h" 46#include "MagickCore/blob.h" 47#include "MagickCore/blob-private.h" 48#include "MagickCore/cache.h" 49#include "MagickCore/cache-private.h" 50#include "MagickCore/cache-view.h" 51#include "MagickCore/client.h" 52#include "MagickCore/color.h" 53#include "MagickCore/color-private.h" 54#include "MagickCore/colorspace.h" 55#include "MagickCore/colorspace-private.h" 56#include "MagickCore/composite.h" 57#include "MagickCore/composite-private.h" 58#include "MagickCore/compress.h" 59#include "MagickCore/constitute.h" 60#include "MagickCore/display.h" 61#include "MagickCore/draw.h" 62#include "MagickCore/enhance.h" 63#include "MagickCore/exception.h" 64#include "MagickCore/exception-private.h" 65#include "MagickCore/feature.h" 66#include "MagickCore/gem.h" 67#include "MagickCore/geometry.h" 68#include "MagickCore/list.h" 69#include "MagickCore/image-private.h" 70#include "MagickCore/magic.h" 71#include "MagickCore/magick.h" 72#include "MagickCore/matrix.h" 73#include "MagickCore/memory_.h" 74#include "MagickCore/module.h" 75#include "MagickCore/monitor.h" 76#include "MagickCore/monitor-private.h" 77#include "MagickCore/morphology-private.h" 78#include "MagickCore/option.h" 79#include "MagickCore/paint.h" 80#include "MagickCore/pixel-accessor.h" 81#include "MagickCore/profile.h" 82#include "MagickCore/quantize.h" 83#include "MagickCore/quantum-private.h" 84#include "MagickCore/random_.h" 85#include "MagickCore/resource_.h" 86#include "MagickCore/segment.h" 87#include "MagickCore/semaphore.h" 88#include "MagickCore/signature-private.h" 89#include "MagickCore/string_.h" 90#include "MagickCore/thread-private.h" 91#include "MagickCore/timer.h" 92#include "MagickCore/utility.h" 93#include "MagickCore/version.h" 94 95/* 96%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 97% % 98% % 99% % 100% C a n n y E d g e I m a g e % 101% % 102% % 103% % 104%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 105% 106% CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of 107% edges in images. 108% 109% The format of the EdgeImage method is: 110% 111% Image *CannyEdgeImage(const Image *image,const double radius, 112% const double sigma,const double lower_percent, 113% const double upper_percent,ExceptionInfo *exception) 114% 115% A description of each parameter follows: 116% 117% o image: the image. 118% 119% o radius: the radius of the gaussian smoothing filter. 120% 121% o sigma: the sigma of the gaussian smoothing filter. 122% 123% o lower_precent: percentage of edge pixels in the lower threshold. 124% 125% o upper_percent: percentage of edge pixels in the upper threshold. 126% 127% o exception: return any errors or warnings in this structure. 128% 129*/ 130 131typedef struct _CannyInfo 132{ 133 double 134 magnitude, 135 intensity; 136 137 int 138 orientation; 139 140 ssize_t 141 x, 142 y; 143} CannyInfo; 144 145static inline MagickBooleanType IsAuthenticPixel(const Image *image, 146 const ssize_t x,const ssize_t y) 147{ 148 if ((x < 0) || (x >= (ssize_t) image->columns)) 149 return(MagickFalse); 150 if ((y < 0) || (y >= (ssize_t) image->rows)) 151 return(MagickFalse); 152 return(MagickTrue); 153} 154 155static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view, 156 MatrixInfo *canny_cache,const ssize_t x,const ssize_t y, 157 const double lower_threshold,ExceptionInfo *exception) 158{ 159 CannyInfo 160 edge, 161 pixel; 162 163 MagickBooleanType 164 status; 165 166 register Quantum 167 *q; 168 169 register ssize_t 170 i; 171 172 q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception); 173 if (q == (Quantum *) NULL) 174 return(MagickFalse); 175 *q=QuantumRange; 176 status=SyncCacheViewAuthenticPixels(edge_view,exception); 177 if (status == MagickFalse) 178 return(MagickFalse);; 179 if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse) 180 return(MagickFalse); 181 edge.x=x; 182 edge.y=y; 183 if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse) 184 return(MagickFalse); 185 for (i=1; i != 0; ) 186 { 187 ssize_t 188 v; 189 190 i--; 191 status=GetMatrixElement(canny_cache,i,0,&edge); 192 if (status == MagickFalse) 193 return(MagickFalse); 194 for (v=(-1); v <= 1; v++) 195 { 196 ssize_t 197 u; 198 199 for (u=(-1); u <= 1; u++) 200 { 201 if ((u == 0) && (v == 0)) 202 continue; 203 if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse) 204 continue; 205 /* 206 Not an edge if gradient value is below the lower threshold. 207 */ 208 q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1, 209 exception); 210 if (q == (Quantum *) NULL) 211 return(MagickFalse); 212 status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel); 213 if (status == MagickFalse) 214 return(MagickFalse); 215 if ((GetPixelIntensity(edge_image,q) == 0.0) && 216 (pixel.intensity >= lower_threshold)) 217 { 218 *q=QuantumRange; 219 status=SyncCacheViewAuthenticPixels(edge_view,exception); 220 if (status == MagickFalse) 221 return(MagickFalse); 222 edge.x+=u; 223 edge.y+=v; 224 status=SetMatrixElement(canny_cache,i,0,&edge); 225 if (status == MagickFalse) 226 return(MagickFalse); 227 i++; 228 } 229 } 230 } 231 } 232 return(MagickTrue); 233} 234 235MagickExport Image *CannyEdgeImage(const Image *image,const double radius, 236 const double sigma,const double lower_percent,const double upper_percent, 237 ExceptionInfo *exception) 238{ 239 CacheView 240 *edge_view; 241 242 CannyInfo 243 pixel; 244 245 char 246 geometry[MaxTextExtent]; 247 248 double 249 lower_threshold, 250 max, 251 min, 252 upper_threshold; 253 254 Image 255 *edge_image; 256 257 KernelInfo 258 *kernel_info; 259 260 MagickBooleanType 261 status; 262 263 MatrixInfo 264 *canny_cache; 265 266 ssize_t 267 y; 268 269 assert(image != (const Image *) NULL); 270 assert(image->signature == MagickSignature); 271 if (image->debug != MagickFalse) 272 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 273 assert(exception != (ExceptionInfo *) NULL); 274 assert(exception->signature == MagickSignature); 275 /* 276 Filter out noise. 277 */ 278 (void) FormatLocaleString(geometry,MaxTextExtent, 279 "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma); 280 kernel_info=AcquireKernelInfo(geometry); 281 if (kernel_info == (KernelInfo *) NULL) 282 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed"); 283 edge_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info, 284 UndefinedCompositeOp,0.0,exception); 285 kernel_info=DestroyKernelInfo(kernel_info); 286 if (edge_image == (Image *) NULL) 287 return((Image *) NULL); 288 if (SetImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse) 289 { 290 edge_image=DestroyImage(edge_image); 291 return((Image *) NULL); 292 } 293 /* 294 Find the intensity gradient of the image. 295 */ 296 canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows, 297 sizeof(CannyInfo),exception); 298 if (canny_cache == (MatrixInfo *) NULL) 299 { 300 edge_image=DestroyImage(edge_image); 301 return((Image *) NULL); 302 } 303 status=MagickTrue; 304 edge_view=AcquireVirtualCacheView(edge_image,exception); 305#if defined(MAGICKCORE_OPENMP_SUPPORT) 306 #pragma omp parallel for schedule(static,4) shared(status) \ 307 magick_threads(edge_image,edge_image,edge_image->rows,1) 308#endif 309 for (y=0; y < (ssize_t) edge_image->rows; y++) 310 { 311 register const Quantum 312 *restrict p; 313 314 register ssize_t 315 x; 316 317 if (status == MagickFalse) 318 continue; 319 p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2, 320 exception); 321 if (p == (const Quantum *) NULL) 322 { 323 status=MagickFalse; 324 continue; 325 } 326 for (x=0; x < (ssize_t) edge_image->columns; x++) 327 { 328 CannyInfo 329 pixel; 330 331 double 332 dx, 333 dy; 334 335 register const Quantum 336 *restrict kernel_pixels; 337 338 ssize_t 339 v; 340 341 static double 342 Gx[2][2] = 343 { 344 { -1.0, +1.0 }, 345 { -1.0, +1.0 } 346 }, 347 Gy[2][2] = 348 { 349 { +1.0, +1.0 }, 350 { -1.0, -1.0 } 351 }; 352 353 (void) ResetMagickMemory(&pixel,0,sizeof(pixel)); 354 dx=0.0; 355 dy=0.0; 356 kernel_pixels=p; 357 for (v=0; v < 2; v++) 358 { 359 ssize_t 360 u; 361 362 for (u=0; u < 2; u++) 363 { 364 double 365 intensity; 366 367 intensity=GetPixelIntensity(edge_image,kernel_pixels+u); 368 dx+=0.5*Gx[v][u]*intensity; 369 dy+=0.5*Gy[v][u]*intensity; 370 } 371 kernel_pixels+=edge_image->columns+1; 372 } 373 pixel.magnitude=hypot(dx,dy); 374 pixel.orientation=0; 375 if (fabs(dx) > MagickEpsilon) 376 { 377 double 378 slope; 379 380 slope=dy/dx; 381 if (slope < 0.0) 382 { 383 if (slope < -2.41421356237) 384 pixel.orientation=0; 385 else 386 if (slope < -0.414213562373) 387 pixel.orientation=1; 388 else 389 pixel.orientation=2; 390 } 391 else 392 { 393 if (slope > 2.41421356237) 394 pixel.orientation=0; 395 else 396 if (slope > 0.414213562373) 397 pixel.orientation=3; 398 else 399 pixel.orientation=2; 400 } 401 } 402 if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse) 403 continue; 404 p+=GetPixelChannels(edge_image); 405 } 406 } 407 edge_view=DestroyCacheView(edge_view); 408 /* 409 Non-maxima suppression, remove pixels that are not considered to be part 410 of an edge. 411 */ 412 (void) GetMatrixElement(canny_cache,0,0,&pixel); 413 max=pixel.intensity; 414 min=pixel.intensity; 415 edge_view=AcquireAuthenticCacheView(edge_image,exception); 416#if defined(MAGICKCORE_OPENMP_SUPPORT) 417 #pragma omp parallel for schedule(static,4) shared(status) \ 418 magick_threads(edge_image,edge_image,edge_image->rows,1) 419#endif 420 for (y=0; y < (ssize_t) edge_image->rows; y++) 421 { 422 register Quantum 423 *restrict q; 424 425 register ssize_t 426 x; 427 428 if (status == MagickFalse) 429 continue; 430 q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1, 431 exception); 432 if (q == (Quantum *) NULL) 433 { 434 status=MagickFalse; 435 continue; 436 } 437 for (x=0; x < (ssize_t) edge_image->columns; x++) 438 { 439 CannyInfo 440 alpha_pixel, 441 beta_pixel, 442 pixel; 443 444 (void) GetMatrixElement(canny_cache,x,y,&pixel); 445 switch (pixel.orientation) 446 { 447 case 0: 448 default: 449 { 450 /* 451 0 degrees, north and south. 452 */ 453 (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel); 454 (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel); 455 break; 456 } 457 case 1: 458 { 459 /* 460 45 degrees, northwest and southeast. 461 */ 462 (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel); 463 (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel); 464 break; 465 } 466 case 2: 467 { 468 /* 469 90 degrees, east and west. 470 */ 471 (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel); 472 (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel); 473 break; 474 } 475 case 3: 476 { 477 /* 478 135 degrees, northeast and southwest. 479 */ 480 (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel); 481 (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel); 482 break; 483 } 484 } 485 pixel.intensity=pixel.magnitude; 486 if ((pixel.magnitude < alpha_pixel.magnitude) || 487 (pixel.magnitude < beta_pixel.magnitude)) 488 pixel.intensity=0; 489 (void) SetMatrixElement(canny_cache,x,y,&pixel); 490#if defined(MAGICKCORE_OPENMP_SUPPORT) 491 #pragma omp critical (MagickCore_CannyEdgeImage) 492#endif 493 { 494 if (pixel.intensity < min) 495 min=pixel.intensity; 496 if (pixel.intensity > max) 497 max=pixel.intensity; 498 } 499 *q=0; 500 q+=GetPixelChannels(edge_image); 501 } 502 if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse) 503 status=MagickFalse; 504 } 505 edge_view=DestroyCacheView(edge_view); 506 /* 507 Estimate hysteresis threshold. 508 */ 509 lower_threshold=lower_percent*(max-min)+min; 510 upper_threshold=upper_percent*(max-min)+min; 511 /* 512 Hysteresis threshold. 513 */ 514 edge_view=AcquireAuthenticCacheView(edge_image,exception); 515 for (y=0; y < (ssize_t) edge_image->rows; y++) 516 { 517 register ssize_t 518 x; 519 520 if (status == MagickFalse) 521 continue; 522 for (x=0; x < (ssize_t) edge_image->columns; x++) 523 { 524 CannyInfo 525 pixel; 526 527 register const Quantum 528 *restrict p; 529 530 /* 531 Edge if pixel gradient higher than upper threshold. 532 */ 533 p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception); 534 if (p == (const Quantum *) NULL) 535 continue; 536 status=GetMatrixElement(canny_cache,x,y,&pixel); 537 if (status == MagickFalse) 538 continue; 539 if ((GetPixelIntensity(edge_image,p) == 0.0) && 540 (pixel.intensity >= upper_threshold)) 541 status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold, 542 exception); 543 } 544 } 545 edge_view=DestroyCacheView(edge_view); 546 /* 547 Free resources. 548 */ 549 canny_cache=DestroyMatrixInfo(canny_cache); 550 return(edge_image); 551} 552 553/* 554%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 555% % 556% % 557% % 558% H o u g h L i n e s I m a g e % 559% % 560% % 561% % 562%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 563% 564% HoughLinesImage() identifies lines in an image. 565% 566% The format of the HoughLinesImage method is: 567% 568% Image *HoughLinesImage(const Image *image,const size_t width, 569% const size_t height,const size_t threshold,ExceptionInfo *exception) 570% 571% A description of each parameter follows: 572% 573% o image: the image. 574% 575% o width, height: find line pairs as local maxima in this neighborhood. 576% 577% o threshold: the line count threshold. 578% 579% o exception: return any errors or warnings in this structure. 580% 581*/ 582MagickExport Image *HoughLinesImage(const Image *image,const size_t width, 583 const size_t height,const size_t threshold,ExceptionInfo *exception) 584{ 585 return((Image *) NULL); 586} 587 588/* 589%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 590% % 591% % 592% % 593% G e t I m a g e F e a t u r e s % 594% % 595% % 596% % 597%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 598% 599% GetImageFeatures() returns features for each channel in the image in 600% each of four directions (horizontal, vertical, left and right diagonals) 601% for the specified distance. The features include the angular second 602% moment, contrast, correlation, sum of squares: variance, inverse difference 603% moment, sum average, sum varience, sum entropy, entropy, difference variance,% difference entropy, information measures of correlation 1, information 604% measures of correlation 2, and maximum correlation coefficient. You can 605% access the red channel contrast, for example, like this: 606% 607% channel_features=GetImageFeatures(image,1,exception); 608% contrast=channel_features[RedPixelChannel].contrast[0]; 609% 610% Use MagickRelinquishMemory() to free the features buffer. 611% 612% The format of the GetImageFeatures method is: 613% 614% ChannelFeatures *GetImageFeatures(const Image *image, 615% const size_t distance,ExceptionInfo *exception) 616% 617% A description of each parameter follows: 618% 619% o image: the image. 620% 621% o distance: the distance. 622% 623% o exception: return any errors or warnings in this structure. 624% 625*/ 626 627static inline ssize_t MagickAbsoluteValue(const ssize_t x) 628{ 629 if (x < 0) 630 return(-x); 631 return(x); 632} 633 634static inline double MagickLog10(const double x) 635{ 636#define Log10Epsilon (1.0e-11) 637 638 if (fabs(x) < Log10Epsilon) 639 return(log10(Log10Epsilon)); 640 return(log10(fabs(x))); 641} 642 643MagickExport ChannelFeatures *GetImageFeatures(const Image *image, 644 const size_t distance,ExceptionInfo *exception) 645{ 646 typedef struct _ChannelStatistics 647 { 648 PixelInfo 649 direction[4]; /* horizontal, vertical, left and right diagonals */ 650 } ChannelStatistics; 651 652 CacheView 653 *image_view; 654 655 ChannelFeatures 656 *channel_features; 657 658 ChannelStatistics 659 **cooccurrence, 660 correlation, 661 *density_x, 662 *density_xy, 663 *density_y, 664 entropy_x, 665 entropy_xy, 666 entropy_xy1, 667 entropy_xy2, 668 entropy_y, 669 mean, 670 **Q, 671 *sum, 672 sum_squares, 673 variance; 674 675 PixelPacket 676 gray, 677 *grays; 678 679 MagickBooleanType 680 status; 681 682 register ssize_t 683 i; 684 685 size_t 686 length; 687 688 ssize_t 689 y; 690 691 unsigned int 692 number_grays; 693 694 assert(image != (Image *) NULL); 695 assert(image->signature == MagickSignature); 696 if (image->debug != MagickFalse) 697 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 698 if ((image->columns < (distance+1)) || (image->rows < (distance+1))) 699 return((ChannelFeatures *) NULL); 700 length=CompositeChannels+1UL; 701 channel_features=(ChannelFeatures *) AcquireQuantumMemory(length, 702 sizeof(*channel_features)); 703 if (channel_features == (ChannelFeatures *) NULL) 704 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed"); 705 (void) ResetMagickMemory(channel_features,0,length* 706 sizeof(*channel_features)); 707 /* 708 Form grays. 709 */ 710 grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays)); 711 if (grays == (PixelPacket *) NULL) 712 { 713 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 714 channel_features); 715 (void) ThrowMagickException(exception,GetMagickModule(), 716 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 717 return(channel_features); 718 } 719 for (i=0; i <= (ssize_t) MaxMap; i++) 720 { 721 grays[i].red=(~0U); 722 grays[i].green=(~0U); 723 grays[i].blue=(~0U); 724 grays[i].alpha=(~0U); 725 grays[i].black=(~0U); 726 } 727 status=MagickTrue; 728 image_view=AcquireVirtualCacheView(image,exception); 729#if defined(MAGICKCORE_OPENMP_SUPPORT) 730 #pragma omp parallel for schedule(static,4) shared(status) \ 731 magick_threads(image,image,image->rows,1) 732#endif 733 for (y=0; y < (ssize_t) image->rows; y++) 734 { 735 register const Quantum 736 *restrict p; 737 738 register ssize_t 739 x; 740 741 if (status == MagickFalse) 742 continue; 743 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 744 if (p == (const Quantum *) NULL) 745 { 746 status=MagickFalse; 747 continue; 748 } 749 for (x=0; x < (ssize_t) image->columns; x++) 750 { 751 grays[ScaleQuantumToMap(GetPixelRed(image,p))].red= 752 ScaleQuantumToMap(GetPixelRed(image,p)); 753 grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green= 754 ScaleQuantumToMap(GetPixelGreen(image,p)); 755 grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue= 756 ScaleQuantumToMap(GetPixelBlue(image,p)); 757 if (image->colorspace == CMYKColorspace) 758 grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black= 759 ScaleQuantumToMap(GetPixelBlack(image,p)); 760 if (image->alpha_trait == BlendPixelTrait) 761 grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha= 762 ScaleQuantumToMap(GetPixelAlpha(image,p)); 763 p+=GetPixelChannels(image); 764 } 765 } 766 image_view=DestroyCacheView(image_view); 767 if (status == MagickFalse) 768 { 769 grays=(PixelPacket *) RelinquishMagickMemory(grays); 770 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 771 channel_features); 772 return(channel_features); 773 } 774 (void) ResetMagickMemory(&gray,0,sizeof(gray)); 775 for (i=0; i <= (ssize_t) MaxMap; i++) 776 { 777 if (grays[i].red != ~0U) 778 grays[gray.red++].red=grays[i].red; 779 if (grays[i].green != ~0U) 780 grays[gray.green++].green=grays[i].green; 781 if (grays[i].blue != ~0U) 782 grays[gray.blue++].blue=grays[i].blue; 783 if (image->colorspace == CMYKColorspace) 784 if (grays[i].black != ~0U) 785 grays[gray.black++].black=grays[i].black; 786 if (image->alpha_trait == BlendPixelTrait) 787 if (grays[i].alpha != ~0U) 788 grays[gray.alpha++].alpha=grays[i].alpha; 789 } 790 /* 791 Allocate spatial dependence matrix. 792 */ 793 number_grays=gray.red; 794 if (gray.green > number_grays) 795 number_grays=gray.green; 796 if (gray.blue > number_grays) 797 number_grays=gray.blue; 798 if (image->colorspace == CMYKColorspace) 799 if (gray.black > number_grays) 800 number_grays=gray.black; 801 if (image->alpha_trait == BlendPixelTrait) 802 if (gray.alpha > number_grays) 803 number_grays=gray.alpha; 804 cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays, 805 sizeof(*cooccurrence)); 806 density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 807 sizeof(*density_x)); 808 density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 809 sizeof(*density_xy)); 810 density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 811 sizeof(*density_y)); 812 Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q)); 813 sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum)); 814 if ((cooccurrence == (ChannelStatistics **) NULL) || 815 (density_x == (ChannelStatistics *) NULL) || 816 (density_xy == (ChannelStatistics *) NULL) || 817 (density_y == (ChannelStatistics *) NULL) || 818 (Q == (ChannelStatistics **) NULL) || 819 (sum == (ChannelStatistics *) NULL)) 820 { 821 if (Q != (ChannelStatistics **) NULL) 822 { 823 for (i=0; i < (ssize_t) number_grays; i++) 824 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 825 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 826 } 827 if (sum != (ChannelStatistics *) NULL) 828 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 829 if (density_y != (ChannelStatistics *) NULL) 830 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 831 if (density_xy != (ChannelStatistics *) NULL) 832 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 833 if (density_x != (ChannelStatistics *) NULL) 834 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 835 if (cooccurrence != (ChannelStatistics **) NULL) 836 { 837 for (i=0; i < (ssize_t) number_grays; i++) 838 cooccurrence[i]=(ChannelStatistics *) 839 RelinquishMagickMemory(cooccurrence[i]); 840 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory( 841 cooccurrence); 842 } 843 grays=(PixelPacket *) RelinquishMagickMemory(grays); 844 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 845 channel_features); 846 (void) ThrowMagickException(exception,GetMagickModule(), 847 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 848 return(channel_features); 849 } 850 (void) ResetMagickMemory(&correlation,0,sizeof(correlation)); 851 (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x)); 852 (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy)); 853 (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y)); 854 (void) ResetMagickMemory(&mean,0,sizeof(mean)); 855 (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum)); 856 (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares)); 857 (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy)); 858 (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x)); 859 (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy)); 860 (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1)); 861 (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2)); 862 (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y)); 863 (void) ResetMagickMemory(&variance,0,sizeof(variance)); 864 for (i=0; i < (ssize_t) number_grays; i++) 865 { 866 cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays, 867 sizeof(**cooccurrence)); 868 Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q)); 869 if ((cooccurrence[i] == (ChannelStatistics *) NULL) || 870 (Q[i] == (ChannelStatistics *) NULL)) 871 break; 872 (void) ResetMagickMemory(cooccurrence[i],0,number_grays* 873 sizeof(**cooccurrence)); 874 (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q)); 875 } 876 if (i < (ssize_t) number_grays) 877 { 878 for (i--; i >= 0; i--) 879 { 880 if (Q[i] != (ChannelStatistics *) NULL) 881 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 882 if (cooccurrence[i] != (ChannelStatistics *) NULL) 883 cooccurrence[i]=(ChannelStatistics *) 884 RelinquishMagickMemory(cooccurrence[i]); 885 } 886 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 887 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 888 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 889 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 890 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 891 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 892 grays=(PixelPacket *) RelinquishMagickMemory(grays); 893 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 894 channel_features); 895 (void) ThrowMagickException(exception,GetMagickModule(), 896 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 897 return(channel_features); 898 } 899 /* 900 Initialize spatial dependence matrix. 901 */ 902 status=MagickTrue; 903 image_view=AcquireVirtualCacheView(image,exception); 904 for (y=0; y < (ssize_t) image->rows; y++) 905 { 906 register const Quantum 907 *restrict p; 908 909 register ssize_t 910 x; 911 912 ssize_t 913 i, 914 offset, 915 u, 916 v; 917 918 if (status == MagickFalse) 919 continue; 920 p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+ 921 2*distance,distance+2,exception); 922 if (p == (const Quantum *) NULL) 923 { 924 status=MagickFalse; 925 continue; 926 } 927 p+=distance*GetPixelChannels(image);; 928 for (x=0; x < (ssize_t) image->columns; x++) 929 { 930 for (i=0; i < 4; i++) 931 { 932 switch (i) 933 { 934 case 0: 935 default: 936 { 937 /* 938 Horizontal adjacency. 939 */ 940 offset=(ssize_t) distance; 941 break; 942 } 943 case 1: 944 { 945 /* 946 Vertical adjacency. 947 */ 948 offset=(ssize_t) (image->columns+2*distance); 949 break; 950 } 951 case 2: 952 { 953 /* 954 Right diagonal adjacency. 955 */ 956 offset=(ssize_t) ((image->columns+2*distance)-distance); 957 break; 958 } 959 case 3: 960 { 961 /* 962 Left diagonal adjacency. 963 */ 964 offset=(ssize_t) ((image->columns+2*distance)+distance); 965 break; 966 } 967 } 968 u=0; 969 v=0; 970 while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p))) 971 u++; 972 while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image)))) 973 v++; 974 cooccurrence[u][v].direction[i].red++; 975 cooccurrence[v][u].direction[i].red++; 976 u=0; 977 v=0; 978 while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p))) 979 u++; 980 while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image)))) 981 v++; 982 cooccurrence[u][v].direction[i].green++; 983 cooccurrence[v][u].direction[i].green++; 984 u=0; 985 v=0; 986 while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p))) 987 u++; 988 while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image)))) 989 v++; 990 cooccurrence[u][v].direction[i].blue++; 991 cooccurrence[v][u].direction[i].blue++; 992 if (image->colorspace == CMYKColorspace) 993 { 994 u=0; 995 v=0; 996 while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p))) 997 u++; 998 while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image)))) 999 v++; 1000 cooccurrence[u][v].direction[i].black++; 1001 cooccurrence[v][u].direction[i].black++; 1002 } 1003 if (image->alpha_trait == BlendPixelTrait) 1004 { 1005 u=0; 1006 v=0; 1007 while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p))) 1008 u++; 1009 while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image)))) 1010 v++; 1011 cooccurrence[u][v].direction[i].alpha++; 1012 cooccurrence[v][u].direction[i].alpha++; 1013 } 1014 } 1015 p+=GetPixelChannels(image); 1016 } 1017 } 1018 grays=(PixelPacket *) RelinquishMagickMemory(grays); 1019 image_view=DestroyCacheView(image_view); 1020 if (status == MagickFalse) 1021 { 1022 for (i=0; i < (ssize_t) number_grays; i++) 1023 cooccurrence[i]=(ChannelStatistics *) 1024 RelinquishMagickMemory(cooccurrence[i]); 1025 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 1026 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 1027 channel_features); 1028 (void) ThrowMagickException(exception,GetMagickModule(), 1029 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 1030 return(channel_features); 1031 } 1032 /* 1033 Normalize spatial dependence matrix. 1034 */ 1035 for (i=0; i < 4; i++) 1036 { 1037 double 1038 normalize; 1039 1040 register ssize_t 1041 y; 1042 1043 switch (i) 1044 { 1045 case 0: 1046 default: 1047 { 1048 /* 1049 Horizontal adjacency. 1050 */ 1051 normalize=2.0*image->rows*(image->columns-distance); 1052 break; 1053 } 1054 case 1: 1055 { 1056 /* 1057 Vertical adjacency. 1058 */ 1059 normalize=2.0*(image->rows-distance)*image->columns; 1060 break; 1061 } 1062 case 2: 1063 { 1064 /* 1065 Right diagonal adjacency. 1066 */ 1067 normalize=2.0*(image->rows-distance)*(image->columns-distance); 1068 break; 1069 } 1070 case 3: 1071 { 1072 /* 1073 Left diagonal adjacency. 1074 */ 1075 normalize=2.0*(image->rows-distance)*(image->columns-distance); 1076 break; 1077 } 1078 } 1079 normalize=PerceptibleReciprocal(normalize); 1080 for (y=0; y < (ssize_t) number_grays; y++) 1081 { 1082 register ssize_t 1083 x; 1084 1085 for (x=0; x < (ssize_t) number_grays; x++) 1086 { 1087 cooccurrence[x][y].direction[i].red*=normalize; 1088 cooccurrence[x][y].direction[i].green*=normalize; 1089 cooccurrence[x][y].direction[i].blue*=normalize; 1090 if (image->colorspace == CMYKColorspace) 1091 cooccurrence[x][y].direction[i].black*=normalize; 1092 if (image->alpha_trait == BlendPixelTrait) 1093 cooccurrence[x][y].direction[i].alpha*=normalize; 1094 } 1095 } 1096 } 1097 /* 1098 Compute texture features. 1099 */ 1100#if defined(MAGICKCORE_OPENMP_SUPPORT) 1101 #pragma omp parallel for schedule(static,4) shared(status) \ 1102 magick_threads(image,image,number_grays,1) 1103#endif 1104 for (i=0; i < 4; i++) 1105 { 1106 register ssize_t 1107 y; 1108 1109 for (y=0; y < (ssize_t) number_grays; y++) 1110 { 1111 register ssize_t 1112 x; 1113 1114 for (x=0; x < (ssize_t) number_grays; x++) 1115 { 1116 /* 1117 Angular second moment: measure of homogeneity of the image. 1118 */ 1119 channel_features[RedPixelChannel].angular_second_moment[i]+= 1120 cooccurrence[x][y].direction[i].red* 1121 cooccurrence[x][y].direction[i].red; 1122 channel_features[GreenPixelChannel].angular_second_moment[i]+= 1123 cooccurrence[x][y].direction[i].green* 1124 cooccurrence[x][y].direction[i].green; 1125 channel_features[BluePixelChannel].angular_second_moment[i]+= 1126 cooccurrence[x][y].direction[i].blue* 1127 cooccurrence[x][y].direction[i].blue; 1128 if (image->colorspace == CMYKColorspace) 1129 channel_features[BlackPixelChannel].angular_second_moment[i]+= 1130 cooccurrence[x][y].direction[i].black* 1131 cooccurrence[x][y].direction[i].black; 1132 if (image->alpha_trait == BlendPixelTrait) 1133 channel_features[AlphaPixelChannel].angular_second_moment[i]+= 1134 cooccurrence[x][y].direction[i].alpha* 1135 cooccurrence[x][y].direction[i].alpha; 1136 /* 1137 Correlation: measure of linear-dependencies in the image. 1138 */ 1139 sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red; 1140 sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green; 1141 sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1142 if (image->colorspace == CMYKColorspace) 1143 sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black; 1144 if (image->alpha_trait == BlendPixelTrait) 1145 sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha; 1146 correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red; 1147 correlation.direction[i].green+=x*y* 1148 cooccurrence[x][y].direction[i].green; 1149 correlation.direction[i].blue+=x*y* 1150 cooccurrence[x][y].direction[i].blue; 1151 if (image->colorspace == CMYKColorspace) 1152 correlation.direction[i].black+=x*y* 1153 cooccurrence[x][y].direction[i].black; 1154 if (image->alpha_trait == BlendPixelTrait) 1155 correlation.direction[i].alpha+=x*y* 1156 cooccurrence[x][y].direction[i].alpha; 1157 /* 1158 Inverse Difference Moment. 1159 */ 1160 channel_features[RedPixelChannel].inverse_difference_moment[i]+= 1161 cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1); 1162 channel_features[GreenPixelChannel].inverse_difference_moment[i]+= 1163 cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1); 1164 channel_features[BluePixelChannel].inverse_difference_moment[i]+= 1165 cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1); 1166 if (image->colorspace == CMYKColorspace) 1167 channel_features[BlackPixelChannel].inverse_difference_moment[i]+= 1168 cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1); 1169 if (image->alpha_trait == BlendPixelTrait) 1170 channel_features[AlphaPixelChannel].inverse_difference_moment[i]+= 1171 cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1); 1172 /* 1173 Sum average. 1174 */ 1175 density_xy[y+x+2].direction[i].red+= 1176 cooccurrence[x][y].direction[i].red; 1177 density_xy[y+x+2].direction[i].green+= 1178 cooccurrence[x][y].direction[i].green; 1179 density_xy[y+x+2].direction[i].blue+= 1180 cooccurrence[x][y].direction[i].blue; 1181 if (image->colorspace == CMYKColorspace) 1182 density_xy[y+x+2].direction[i].black+= 1183 cooccurrence[x][y].direction[i].black; 1184 if (image->alpha_trait == BlendPixelTrait) 1185 density_xy[y+x+2].direction[i].alpha+= 1186 cooccurrence[x][y].direction[i].alpha; 1187 /* 1188 Entropy. 1189 */ 1190 channel_features[RedPixelChannel].entropy[i]-= 1191 cooccurrence[x][y].direction[i].red* 1192 MagickLog10(cooccurrence[x][y].direction[i].red); 1193 channel_features[GreenPixelChannel].entropy[i]-= 1194 cooccurrence[x][y].direction[i].green* 1195 MagickLog10(cooccurrence[x][y].direction[i].green); 1196 channel_features[BluePixelChannel].entropy[i]-= 1197 cooccurrence[x][y].direction[i].blue* 1198 MagickLog10(cooccurrence[x][y].direction[i].blue); 1199 if (image->colorspace == CMYKColorspace) 1200 channel_features[BlackPixelChannel].entropy[i]-= 1201 cooccurrence[x][y].direction[i].black* 1202 MagickLog10(cooccurrence[x][y].direction[i].black); 1203 if (image->alpha_trait == BlendPixelTrait) 1204 channel_features[AlphaPixelChannel].entropy[i]-= 1205 cooccurrence[x][y].direction[i].alpha* 1206 MagickLog10(cooccurrence[x][y].direction[i].alpha); 1207 /* 1208 Information Measures of Correlation. 1209 */ 1210 density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red; 1211 density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green; 1212 density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1213 if (image->alpha_trait == BlendPixelTrait) 1214 density_x[x].direction[i].alpha+= 1215 cooccurrence[x][y].direction[i].alpha; 1216 if (image->colorspace == CMYKColorspace) 1217 density_x[x].direction[i].black+= 1218 cooccurrence[x][y].direction[i].black; 1219 density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red; 1220 density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green; 1221 density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1222 if (image->colorspace == CMYKColorspace) 1223 density_y[y].direction[i].black+= 1224 cooccurrence[x][y].direction[i].black; 1225 if (image->alpha_trait == BlendPixelTrait) 1226 density_y[y].direction[i].alpha+= 1227 cooccurrence[x][y].direction[i].alpha; 1228 } 1229 mean.direction[i].red+=y*sum[y].direction[i].red; 1230 sum_squares.direction[i].red+=y*y*sum[y].direction[i].red; 1231 mean.direction[i].green+=y*sum[y].direction[i].green; 1232 sum_squares.direction[i].green+=y*y*sum[y].direction[i].green; 1233 mean.direction[i].blue+=y*sum[y].direction[i].blue; 1234 sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue; 1235 if (image->colorspace == CMYKColorspace) 1236 { 1237 mean.direction[i].black+=y*sum[y].direction[i].black; 1238 sum_squares.direction[i].black+=y*y*sum[y].direction[i].black; 1239 } 1240 if (image->alpha_trait == BlendPixelTrait) 1241 { 1242 mean.direction[i].alpha+=y*sum[y].direction[i].alpha; 1243 sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha; 1244 } 1245 } 1246 /* 1247 Correlation: measure of linear-dependencies in the image. 1248 */ 1249 channel_features[RedPixelChannel].correlation[i]= 1250 (correlation.direction[i].red-mean.direction[i].red* 1251 mean.direction[i].red)/(sqrt(sum_squares.direction[i].red- 1252 (mean.direction[i].red*mean.direction[i].red))*sqrt( 1253 sum_squares.direction[i].red-(mean.direction[i].red* 1254 mean.direction[i].red))); 1255 channel_features[GreenPixelChannel].correlation[i]= 1256 (correlation.direction[i].green-mean.direction[i].green* 1257 mean.direction[i].green)/(sqrt(sum_squares.direction[i].green- 1258 (mean.direction[i].green*mean.direction[i].green))*sqrt( 1259 sum_squares.direction[i].green-(mean.direction[i].green* 1260 mean.direction[i].green))); 1261 channel_features[BluePixelChannel].correlation[i]= 1262 (correlation.direction[i].blue-mean.direction[i].blue* 1263 mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue- 1264 (mean.direction[i].blue*mean.direction[i].blue))*sqrt( 1265 sum_squares.direction[i].blue-(mean.direction[i].blue* 1266 mean.direction[i].blue))); 1267 if (image->colorspace == CMYKColorspace) 1268 channel_features[BlackPixelChannel].correlation[i]= 1269 (correlation.direction[i].black-mean.direction[i].black* 1270 mean.direction[i].black)/(sqrt(sum_squares.direction[i].black- 1271 (mean.direction[i].black*mean.direction[i].black))*sqrt( 1272 sum_squares.direction[i].black-(mean.direction[i].black* 1273 mean.direction[i].black))); 1274 if (image->alpha_trait == BlendPixelTrait) 1275 channel_features[AlphaPixelChannel].correlation[i]= 1276 (correlation.direction[i].alpha-mean.direction[i].alpha* 1277 mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha- 1278 (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt( 1279 sum_squares.direction[i].alpha-(mean.direction[i].alpha* 1280 mean.direction[i].alpha))); 1281 } 1282 /* 1283 Compute more texture features. 1284 */ 1285#if defined(MAGICKCORE_OPENMP_SUPPORT) 1286 #pragma omp parallel for schedule(static,4) shared(status) \ 1287 magick_threads(image,image,number_grays,1) 1288#endif 1289 for (i=0; i < 4; i++) 1290 { 1291 register ssize_t 1292 x; 1293 1294 for (x=2; x < (ssize_t) (2*number_grays); x++) 1295 { 1296 /* 1297 Sum average. 1298 */ 1299 channel_features[RedPixelChannel].sum_average[i]+= 1300 x*density_xy[x].direction[i].red; 1301 channel_features[GreenPixelChannel].sum_average[i]+= 1302 x*density_xy[x].direction[i].green; 1303 channel_features[BluePixelChannel].sum_average[i]+= 1304 x*density_xy[x].direction[i].blue; 1305 if (image->colorspace == CMYKColorspace) 1306 channel_features[BlackPixelChannel].sum_average[i]+= 1307 x*density_xy[x].direction[i].black; 1308 if (image->alpha_trait == BlendPixelTrait) 1309 channel_features[AlphaPixelChannel].sum_average[i]+= 1310 x*density_xy[x].direction[i].alpha; 1311 /* 1312 Sum entropy. 1313 */ 1314 channel_features[RedPixelChannel].sum_entropy[i]-= 1315 density_xy[x].direction[i].red* 1316 MagickLog10(density_xy[x].direction[i].red); 1317 channel_features[GreenPixelChannel].sum_entropy[i]-= 1318 density_xy[x].direction[i].green* 1319 MagickLog10(density_xy[x].direction[i].green); 1320 channel_features[BluePixelChannel].sum_entropy[i]-= 1321 density_xy[x].direction[i].blue* 1322 MagickLog10(density_xy[x].direction[i].blue); 1323 if (image->colorspace == CMYKColorspace) 1324 channel_features[BlackPixelChannel].sum_entropy[i]-= 1325 density_xy[x].direction[i].black* 1326 MagickLog10(density_xy[x].direction[i].black); 1327 if (image->alpha_trait == BlendPixelTrait) 1328 channel_features[AlphaPixelChannel].sum_entropy[i]-= 1329 density_xy[x].direction[i].alpha* 1330 MagickLog10(density_xy[x].direction[i].alpha); 1331 /* 1332 Sum variance. 1333 */ 1334 channel_features[RedPixelChannel].sum_variance[i]+= 1335 (x-channel_features[RedPixelChannel].sum_entropy[i])* 1336 (x-channel_features[RedPixelChannel].sum_entropy[i])* 1337 density_xy[x].direction[i].red; 1338 channel_features[GreenPixelChannel].sum_variance[i]+= 1339 (x-channel_features[GreenPixelChannel].sum_entropy[i])* 1340 (x-channel_features[GreenPixelChannel].sum_entropy[i])* 1341 density_xy[x].direction[i].green; 1342 channel_features[BluePixelChannel].sum_variance[i]+= 1343 (x-channel_features[BluePixelChannel].sum_entropy[i])* 1344 (x-channel_features[BluePixelChannel].sum_entropy[i])* 1345 density_xy[x].direction[i].blue; 1346 if (image->colorspace == CMYKColorspace) 1347 channel_features[BlackPixelChannel].sum_variance[i]+= 1348 (x-channel_features[BlackPixelChannel].sum_entropy[i])* 1349 (x-channel_features[BlackPixelChannel].sum_entropy[i])* 1350 density_xy[x].direction[i].black; 1351 if (image->alpha_trait == BlendPixelTrait) 1352 channel_features[AlphaPixelChannel].sum_variance[i]+= 1353 (x-channel_features[AlphaPixelChannel].sum_entropy[i])* 1354 (x-channel_features[AlphaPixelChannel].sum_entropy[i])* 1355 density_xy[x].direction[i].alpha; 1356 } 1357 } 1358 /* 1359 Compute more texture features. 1360 */ 1361#if defined(MAGICKCORE_OPENMP_SUPPORT) 1362 #pragma omp parallel for schedule(static,4) shared(status) \ 1363 magick_threads(image,image,number_grays,1) 1364#endif 1365 for (i=0; i < 4; i++) 1366 { 1367 register ssize_t 1368 y; 1369 1370 for (y=0; y < (ssize_t) number_grays; y++) 1371 { 1372 register ssize_t 1373 x; 1374 1375 for (x=0; x < (ssize_t) number_grays; x++) 1376 { 1377 /* 1378 Sum of Squares: Variance 1379 */ 1380 variance.direction[i].red+=(y-mean.direction[i].red+1)* 1381 (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red; 1382 variance.direction[i].green+=(y-mean.direction[i].green+1)* 1383 (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green; 1384 variance.direction[i].blue+=(y-mean.direction[i].blue+1)* 1385 (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue; 1386 if (image->colorspace == CMYKColorspace) 1387 variance.direction[i].black+=(y-mean.direction[i].black+1)* 1388 (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black; 1389 if (image->alpha_trait == BlendPixelTrait) 1390 variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)* 1391 (y-mean.direction[i].alpha+1)* 1392 cooccurrence[x][y].direction[i].alpha; 1393 /* 1394 Sum average / Difference Variance. 1395 */ 1396 density_xy[MagickAbsoluteValue(y-x)].direction[i].red+= 1397 cooccurrence[x][y].direction[i].red; 1398 density_xy[MagickAbsoluteValue(y-x)].direction[i].green+= 1399 cooccurrence[x][y].direction[i].green; 1400 density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+= 1401 cooccurrence[x][y].direction[i].blue; 1402 if (image->colorspace == CMYKColorspace) 1403 density_xy[MagickAbsoluteValue(y-x)].direction[i].black+= 1404 cooccurrence[x][y].direction[i].black; 1405 if (image->alpha_trait == BlendPixelTrait) 1406 density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+= 1407 cooccurrence[x][y].direction[i].alpha; 1408 /* 1409 Information Measures of Correlation. 1410 */ 1411 entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red* 1412 MagickLog10(cooccurrence[x][y].direction[i].red); 1413 entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green* 1414 MagickLog10(cooccurrence[x][y].direction[i].green); 1415 entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue* 1416 MagickLog10(cooccurrence[x][y].direction[i].blue); 1417 if (image->colorspace == CMYKColorspace) 1418 entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black* 1419 MagickLog10(cooccurrence[x][y].direction[i].black); 1420 if (image->alpha_trait == BlendPixelTrait) 1421 entropy_xy.direction[i].alpha-= 1422 cooccurrence[x][y].direction[i].alpha*MagickLog10( 1423 cooccurrence[x][y].direction[i].alpha); 1424 entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red* 1425 MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red)); 1426 entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green* 1427 MagickLog10(density_x[x].direction[i].green* 1428 density_y[y].direction[i].green)); 1429 entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue* 1430 MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue)); 1431 if (image->colorspace == CMYKColorspace) 1432 entropy_xy1.direction[i].black-=( 1433 cooccurrence[x][y].direction[i].black*MagickLog10( 1434 density_x[x].direction[i].black*density_y[y].direction[i].black)); 1435 if (image->alpha_trait == BlendPixelTrait) 1436 entropy_xy1.direction[i].alpha-=( 1437 cooccurrence[x][y].direction[i].alpha*MagickLog10( 1438 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha)); 1439 entropy_xy2.direction[i].red-=(density_x[x].direction[i].red* 1440 density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red* 1441 density_y[y].direction[i].red)); 1442 entropy_xy2.direction[i].green-=(density_x[x].direction[i].green* 1443 density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green* 1444 density_y[y].direction[i].green)); 1445 entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue* 1446 density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue* 1447 density_y[y].direction[i].blue)); 1448 if (image->colorspace == CMYKColorspace) 1449 entropy_xy2.direction[i].black-=(density_x[x].direction[i].black* 1450 density_y[y].direction[i].black*MagickLog10( 1451 density_x[x].direction[i].black*density_y[y].direction[i].black)); 1452 if (image->alpha_trait == BlendPixelTrait) 1453 entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha* 1454 density_y[y].direction[i].alpha*MagickLog10( 1455 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha)); 1456 } 1457 } 1458 channel_features[RedPixelChannel].variance_sum_of_squares[i]= 1459 variance.direction[i].red; 1460 channel_features[GreenPixelChannel].variance_sum_of_squares[i]= 1461 variance.direction[i].green; 1462 channel_features[BluePixelChannel].variance_sum_of_squares[i]= 1463 variance.direction[i].blue; 1464 if (image->colorspace == CMYKColorspace) 1465 channel_features[BlackPixelChannel].variance_sum_of_squares[i]= 1466 variance.direction[i].black; 1467 if (image->alpha_trait == BlendPixelTrait) 1468 channel_features[AlphaPixelChannel].variance_sum_of_squares[i]= 1469 variance.direction[i].alpha; 1470 } 1471 /* 1472 Compute more texture features. 1473 */ 1474 (void) ResetMagickMemory(&variance,0,sizeof(variance)); 1475 (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares)); 1476#if defined(MAGICKCORE_OPENMP_SUPPORT) 1477 #pragma omp parallel for schedule(static,4) shared(status) \ 1478 magick_threads(image,image,number_grays,1) 1479#endif 1480 for (i=0; i < 4; i++) 1481 { 1482 register ssize_t 1483 x; 1484 1485 for (x=0; x < (ssize_t) number_grays; x++) 1486 { 1487 /* 1488 Difference variance. 1489 */ 1490 variance.direction[i].red+=density_xy[x].direction[i].red; 1491 variance.direction[i].green+=density_xy[x].direction[i].green; 1492 variance.direction[i].blue+=density_xy[x].direction[i].blue; 1493 if (image->colorspace == CMYKColorspace) 1494 variance.direction[i].black+=density_xy[x].direction[i].black; 1495 if (image->alpha_trait == BlendPixelTrait) 1496 variance.direction[i].alpha+=density_xy[x].direction[i].alpha; 1497 sum_squares.direction[i].red+=density_xy[x].direction[i].red* 1498 density_xy[x].direction[i].red; 1499 sum_squares.direction[i].green+=density_xy[x].direction[i].green* 1500 density_xy[x].direction[i].green; 1501 sum_squares.direction[i].blue+=density_xy[x].direction[i].blue* 1502 density_xy[x].direction[i].blue; 1503 if (image->colorspace == CMYKColorspace) 1504 sum_squares.direction[i].black+=density_xy[x].direction[i].black* 1505 density_xy[x].direction[i].black; 1506 if (image->alpha_trait == BlendPixelTrait) 1507 sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha* 1508 density_xy[x].direction[i].alpha; 1509 /* 1510 Difference entropy. 1511 */ 1512 channel_features[RedPixelChannel].difference_entropy[i]-= 1513 density_xy[x].direction[i].red* 1514 MagickLog10(density_xy[x].direction[i].red); 1515 channel_features[GreenPixelChannel].difference_entropy[i]-= 1516 density_xy[x].direction[i].green* 1517 MagickLog10(density_xy[x].direction[i].green); 1518 channel_features[BluePixelChannel].difference_entropy[i]-= 1519 density_xy[x].direction[i].blue* 1520 MagickLog10(density_xy[x].direction[i].blue); 1521 if (image->colorspace == CMYKColorspace) 1522 channel_features[BlackPixelChannel].difference_entropy[i]-= 1523 density_xy[x].direction[i].black* 1524 MagickLog10(density_xy[x].direction[i].black); 1525 if (image->alpha_trait == BlendPixelTrait) 1526 channel_features[AlphaPixelChannel].difference_entropy[i]-= 1527 density_xy[x].direction[i].alpha* 1528 MagickLog10(density_xy[x].direction[i].alpha); 1529 /* 1530 Information Measures of Correlation. 1531 */ 1532 entropy_x.direction[i].red-=(density_x[x].direction[i].red* 1533 MagickLog10(density_x[x].direction[i].red)); 1534 entropy_x.direction[i].green-=(density_x[x].direction[i].green* 1535 MagickLog10(density_x[x].direction[i].green)); 1536 entropy_x.direction[i].blue-=(density_x[x].direction[i].blue* 1537 MagickLog10(density_x[x].direction[i].blue)); 1538 if (image->colorspace == CMYKColorspace) 1539 entropy_x.direction[i].black-=(density_x[x].direction[i].black* 1540 MagickLog10(density_x[x].direction[i].black)); 1541 if (image->alpha_trait == BlendPixelTrait) 1542 entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha* 1543 MagickLog10(density_x[x].direction[i].alpha)); 1544 entropy_y.direction[i].red-=(density_y[x].direction[i].red* 1545 MagickLog10(density_y[x].direction[i].red)); 1546 entropy_y.direction[i].green-=(density_y[x].direction[i].green* 1547 MagickLog10(density_y[x].direction[i].green)); 1548 entropy_y.direction[i].blue-=(density_y[x].direction[i].blue* 1549 MagickLog10(density_y[x].direction[i].blue)); 1550 if (image->colorspace == CMYKColorspace) 1551 entropy_y.direction[i].black-=(density_y[x].direction[i].black* 1552 MagickLog10(density_y[x].direction[i].black)); 1553 if (image->alpha_trait == BlendPixelTrait) 1554 entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha* 1555 MagickLog10(density_y[x].direction[i].alpha)); 1556 } 1557 /* 1558 Difference variance. 1559 */ 1560 channel_features[RedPixelChannel].difference_variance[i]= 1561 (((double) number_grays*number_grays*sum_squares.direction[i].red)- 1562 (variance.direction[i].red*variance.direction[i].red))/ 1563 ((double) number_grays*number_grays*number_grays*number_grays); 1564 channel_features[GreenPixelChannel].difference_variance[i]= 1565 (((double) number_grays*number_grays*sum_squares.direction[i].green)- 1566 (variance.direction[i].green*variance.direction[i].green))/ 1567 ((double) number_grays*number_grays*number_grays*number_grays); 1568 channel_features[BluePixelChannel].difference_variance[i]= 1569 (((double) number_grays*number_grays*sum_squares.direction[i].blue)- 1570 (variance.direction[i].blue*variance.direction[i].blue))/ 1571 ((double) number_grays*number_grays*number_grays*number_grays); 1572 if (image->colorspace == CMYKColorspace) 1573 channel_features[BlackPixelChannel].difference_variance[i]= 1574 (((double) number_grays*number_grays*sum_squares.direction[i].black)- 1575 (variance.direction[i].black*variance.direction[i].black))/ 1576 ((double) number_grays*number_grays*number_grays*number_grays); 1577 if (image->alpha_trait == BlendPixelTrait) 1578 channel_features[AlphaPixelChannel].difference_variance[i]= 1579 (((double) number_grays*number_grays*sum_squares.direction[i].alpha)- 1580 (variance.direction[i].alpha*variance.direction[i].alpha))/ 1581 ((double) number_grays*number_grays*number_grays*number_grays); 1582 /* 1583 Information Measures of Correlation. 1584 */ 1585 channel_features[RedPixelChannel].measure_of_correlation_1[i]= 1586 (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/ 1587 (entropy_x.direction[i].red > entropy_y.direction[i].red ? 1588 entropy_x.direction[i].red : entropy_y.direction[i].red); 1589 channel_features[GreenPixelChannel].measure_of_correlation_1[i]= 1590 (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/ 1591 (entropy_x.direction[i].green > entropy_y.direction[i].green ? 1592 entropy_x.direction[i].green : entropy_y.direction[i].green); 1593 channel_features[BluePixelChannel].measure_of_correlation_1[i]= 1594 (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/ 1595 (entropy_x.direction[i].blue > entropy_y.direction[i].blue ? 1596 entropy_x.direction[i].blue : entropy_y.direction[i].blue); 1597 if (image->colorspace == CMYKColorspace) 1598 channel_features[BlackPixelChannel].measure_of_correlation_1[i]= 1599 (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/ 1600 (entropy_x.direction[i].black > entropy_y.direction[i].black ? 1601 entropy_x.direction[i].black : entropy_y.direction[i].black); 1602 if (image->alpha_trait == BlendPixelTrait) 1603 channel_features[AlphaPixelChannel].measure_of_correlation_1[i]= 1604 (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/ 1605 (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ? 1606 entropy_x.direction[i].alpha : entropy_y.direction[i].alpha); 1607 channel_features[RedPixelChannel].measure_of_correlation_2[i]= 1608 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red- 1609 entropy_xy.direction[i].red))))); 1610 channel_features[GreenPixelChannel].measure_of_correlation_2[i]= 1611 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green- 1612 entropy_xy.direction[i].green))))); 1613 channel_features[BluePixelChannel].measure_of_correlation_2[i]= 1614 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue- 1615 entropy_xy.direction[i].blue))))); 1616 if (image->colorspace == CMYKColorspace) 1617 channel_features[BlackPixelChannel].measure_of_correlation_2[i]= 1618 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black- 1619 entropy_xy.direction[i].black))))); 1620 if (image->alpha_trait == BlendPixelTrait) 1621 channel_features[AlphaPixelChannel].measure_of_correlation_2[i]= 1622 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha- 1623 entropy_xy.direction[i].alpha))))); 1624 } 1625 /* 1626 Compute more texture features. 1627 */ 1628#if defined(MAGICKCORE_OPENMP_SUPPORT) 1629 #pragma omp parallel for schedule(static,4) shared(status) \ 1630 magick_threads(image,image,number_grays,1) 1631#endif 1632 for (i=0; i < 4; i++) 1633 { 1634 ssize_t 1635 z; 1636 1637 for (z=0; z < (ssize_t) number_grays; z++) 1638 { 1639 register ssize_t 1640 y; 1641 1642 ChannelStatistics 1643 pixel; 1644 1645 (void) ResetMagickMemory(&pixel,0,sizeof(pixel)); 1646 for (y=0; y < (ssize_t) number_grays; y++) 1647 { 1648 register ssize_t 1649 x; 1650 1651 for (x=0; x < (ssize_t) number_grays; x++) 1652 { 1653 /* 1654 Contrast: amount of local variations present in an image. 1655 */ 1656 if (((y-x) == z) || ((x-y) == z)) 1657 { 1658 pixel.direction[i].red+=cooccurrence[x][y].direction[i].red; 1659 pixel.direction[i].green+=cooccurrence[x][y].direction[i].green; 1660 pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1661 if (image->colorspace == CMYKColorspace) 1662 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black; 1663 if (image->alpha_trait == BlendPixelTrait) 1664 pixel.direction[i].alpha+= 1665 cooccurrence[x][y].direction[i].alpha; 1666 } 1667 /* 1668 Maximum Correlation Coefficient. 1669 */ 1670 Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red* 1671 cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/ 1672 density_y[x].direction[i].red; 1673 Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green* 1674 cooccurrence[y][x].direction[i].green/ 1675 density_x[z].direction[i].green/density_y[x].direction[i].red; 1676 Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue* 1677 cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/ 1678 density_y[x].direction[i].blue; 1679 if (image->colorspace == CMYKColorspace) 1680 Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black* 1681 cooccurrence[y][x].direction[i].black/ 1682 density_x[z].direction[i].black/density_y[x].direction[i].black; 1683 if (image->alpha_trait == BlendPixelTrait) 1684 Q[z][y].direction[i].alpha+= 1685 cooccurrence[z][x].direction[i].alpha* 1686 cooccurrence[y][x].direction[i].alpha/ 1687 density_x[z].direction[i].alpha/ 1688 density_y[x].direction[i].alpha; 1689 } 1690 } 1691 channel_features[RedPixelChannel].contrast[i]+=z*z* 1692 pixel.direction[i].red; 1693 channel_features[GreenPixelChannel].contrast[i]+=z*z* 1694 pixel.direction[i].green; 1695 channel_features[BluePixelChannel].contrast[i]+=z*z* 1696 pixel.direction[i].blue; 1697 if (image->colorspace == CMYKColorspace) 1698 channel_features[BlackPixelChannel].contrast[i]+=z*z* 1699 pixel.direction[i].black; 1700 if (image->alpha_trait == BlendPixelTrait) 1701 channel_features[AlphaPixelChannel].contrast[i]+=z*z* 1702 pixel.direction[i].alpha; 1703 } 1704 /* 1705 Maximum Correlation Coefficient. 1706 Future: return second largest eigenvalue of Q. 1707 */ 1708 channel_features[RedPixelChannel].maximum_correlation_coefficient[i]= 1709 sqrt((double) -1.0); 1710 channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]= 1711 sqrt((double) -1.0); 1712 channel_features[BluePixelChannel].maximum_correlation_coefficient[i]= 1713 sqrt((double) -1.0); 1714 if (image->colorspace == CMYKColorspace) 1715 channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]= 1716 sqrt((double) -1.0); 1717 if (image->alpha_trait == BlendPixelTrait) 1718 channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]= 1719 sqrt((double) -1.0); 1720 } 1721 /* 1722 Relinquish resources. 1723 */ 1724 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 1725 for (i=0; i < (ssize_t) number_grays; i++) 1726 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 1727 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 1728 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 1729 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 1730 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 1731 for (i=0; i < (ssize_t) number_grays; i++) 1732 cooccurrence[i]=(ChannelStatistics *) 1733 RelinquishMagickMemory(cooccurrence[i]); 1734 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 1735 return(channel_features); 1736} 1737