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-2016 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 /* 42 Include declarations. 43 */ 44 #include "MagickCore/studio.h" 45 #include "MagickCore/animate.h" 46 #include "MagickCore/artifact.h" 47 #include "MagickCore/blob.h" 48 #include "MagickCore/blob-private.h" 49 #include "MagickCore/cache.h" 50 #include "MagickCore/cache-private.h" 51 #include "MagickCore/cache-view.h" 52 #include "MagickCore/channel.h" 53 #include "MagickCore/client.h" 54 #include "MagickCore/color.h" 55 #include "MagickCore/color-private.h" 56 #include "MagickCore/colorspace.h" 57 #include "MagickCore/colorspace-private.h" 58 #include "MagickCore/composite.h" 59 #include "MagickCore/composite-private.h" 60 #include "MagickCore/compress.h" 61 #include "MagickCore/constitute.h" 62 #include "MagickCore/display.h" 63 #include "MagickCore/draw.h" 64 #include "MagickCore/enhance.h" 65 #include "MagickCore/exception.h" 66 #include "MagickCore/exception-private.h" 67 #include "MagickCore/feature.h" 68 #include "MagickCore/gem.h" 69 #include "MagickCore/geometry.h" 70 #include "MagickCore/list.h" 71 #include "MagickCore/image-private.h" 72 #include "MagickCore/magic.h" 73 #include "MagickCore/magick.h" 74 #include "MagickCore/matrix.h" 75 #include "MagickCore/memory_.h" 76 #include "MagickCore/module.h" 77 #include "MagickCore/monitor.h" 78 #include "MagickCore/monitor-private.h" 79 #include "MagickCore/morphology-private.h" 80 #include "MagickCore/option.h" 81 #include "MagickCore/paint.h" 82 #include "MagickCore/pixel-accessor.h" 83 #include "MagickCore/profile.h" 84 #include "MagickCore/property.h" 85 #include "MagickCore/quantize.h" 86 #include "MagickCore/quantum-private.h" 87 #include "MagickCore/random_.h" 88 #include "MagickCore/resource_.h" 89 #include "MagickCore/segment.h" 90 #include "MagickCore/semaphore.h" 91 #include "MagickCore/signature-private.h" 92 #include "MagickCore/string_.h" 93 #include "MagickCore/thread-private.h" 94 #include "MagickCore/timer.h" 95 #include "MagickCore/utility.h" 96 #include "MagickCore/version.h" 97 98 /* 100 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 101 % % 102 % % 103 % % 104 % C a n n y E d g e I m a g e % 105 % % 106 % % 107 % % 108 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 109 % 110 % CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of 111 % edges in images. 112 % 113 % The format of the CannyEdgeImage method is: 114 % 115 % Image *CannyEdgeImage(const Image *image,const double radius, 116 % const double sigma,const double lower_percent, 117 % const double upper_percent,ExceptionInfo *exception) 118 % 119 % A description of each parameter follows: 120 % 121 % o image: the image. 122 % 123 % o radius: the radius of the gaussian smoothing filter. 124 % 125 % o sigma: the sigma of the gaussian smoothing filter. 126 % 127 % o lower_precent: percentage of edge pixels in the lower threshold. 128 % 129 % o upper_percent: percentage of edge pixels in the upper threshold. 130 % 131 % o exception: return any errors or warnings in this structure. 132 % 133 */ 134 135 typedef struct _CannyInfo 136 { 137 double 138 magnitude, 139 intensity; 140 141 int 142 orientation; 143 144 ssize_t 145 x, 146 y; 147 } CannyInfo; 148 149 static inline MagickBooleanType IsAuthenticPixel(const Image *image, 150 const ssize_t x,const ssize_t y) 151 { 152 if ((x < 0) || (x >= (ssize_t) image->columns)) 153 return(MagickFalse); 154 if ((y < 0) || (y >= (ssize_t) image->rows)) 155 return(MagickFalse); 156 return(MagickTrue); 157 } 158 159 static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view, 160 MatrixInfo *canny_cache,const ssize_t x,const ssize_t y, 161 const double lower_threshold,ExceptionInfo *exception) 162 { 163 CannyInfo 164 edge, 165 pixel; 166 167 MagickBooleanType 168 status; 169 170 register Quantum 171 *q; 172 173 register ssize_t 174 i; 175 176 q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception); 177 if (q == (Quantum *) NULL) 178 return(MagickFalse); 179 *q=QuantumRange; 180 status=SyncCacheViewAuthenticPixels(edge_view,exception); 181 if (status == MagickFalse) 182 return(MagickFalse);; 183 if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse) 184 return(MagickFalse); 185 edge.x=x; 186 edge.y=y; 187 if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse) 188 return(MagickFalse); 189 for (i=1; i != 0; ) 190 { 191 ssize_t 192 v; 193 194 i--; 195 status=GetMatrixElement(canny_cache,i,0,&edge); 196 if (status == MagickFalse) 197 return(MagickFalse); 198 for (v=(-1); v <= 1; v++) 199 { 200 ssize_t 201 u; 202 203 for (u=(-1); u <= 1; u++) 204 { 205 if ((u == 0) && (v == 0)) 206 continue; 207 if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse) 208 continue; 209 /* 210 Not an edge if gradient value is below the lower threshold. 211 */ 212 q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1, 213 exception); 214 if (q == (Quantum *) NULL) 215 return(MagickFalse); 216 status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel); 217 if (status == MagickFalse) 218 return(MagickFalse); 219 if ((GetPixelIntensity(edge_image,q) == 0.0) && 220 (pixel.intensity >= lower_threshold)) 221 { 222 *q=QuantumRange; 223 status=SyncCacheViewAuthenticPixels(edge_view,exception); 224 if (status == MagickFalse) 225 return(MagickFalse); 226 edge.x+=u; 227 edge.y+=v; 228 status=SetMatrixElement(canny_cache,i,0,&edge); 229 if (status == MagickFalse) 230 return(MagickFalse); 231 i++; 232 } 233 } 234 } 235 } 236 return(MagickTrue); 237 } 238 239 MagickExport Image *CannyEdgeImage(const Image *image,const double radius, 240 const double sigma,const double lower_percent,const double upper_percent, 241 ExceptionInfo *exception) 242 { 243 #define CannyEdgeImageTag "CannyEdge/Image" 244 245 CacheView 246 *edge_view; 247 248 CannyInfo 249 element; 250 251 char 252 geometry[MagickPathExtent]; 253 254 double 255 lower_threshold, 256 max, 257 min, 258 upper_threshold; 259 260 Image 261 *edge_image; 262 263 KernelInfo 264 *kernel_info; 265 266 MagickBooleanType 267 status; 268 269 MagickOffsetType 270 progress; 271 272 MatrixInfo 273 *canny_cache; 274 275 ssize_t 276 y; 277 278 assert(image != (const Image *) NULL); 279 assert(image->signature == MagickCoreSignature); 280 if (image->debug != MagickFalse) 281 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 282 assert(exception != (ExceptionInfo *) NULL); 283 assert(exception->signature == MagickCoreSignature); 284 /* 285 Filter out noise. 286 */ 287 (void) FormatLocaleString(geometry,MagickPathExtent, 288 "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma); 289 kernel_info=AcquireKernelInfo(geometry,exception); 290 if (kernel_info == (KernelInfo *) NULL) 291 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed"); 292 edge_image=ConvolveImage(image, kernel_info, exception); 293 kernel_info=DestroyKernelInfo(kernel_info); 294 if (edge_image == (Image *) NULL) 295 return((Image *) NULL); 296 if (SetImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse) 297 { 298 edge_image=DestroyImage(edge_image); 299 return((Image *) NULL); 300 } 301 (void) SetImageAlphaChannel(edge_image,OffAlphaChannel,exception); 302 /* 303 Find the intensity gradient of the image. 304 */ 305 canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows, 306 sizeof(CannyInfo),exception); 307 if (canny_cache == (MatrixInfo *) NULL) 308 { 309 edge_image=DestroyImage(edge_image); 310 return((Image *) NULL); 311 } 312 status=MagickTrue; 313 edge_view=AcquireVirtualCacheView(edge_image,exception); 314 #if defined(MAGICKCORE_OPENMP_SUPPORT) 315 #pragma omp parallel for schedule(static,4) shared(status) \ 316 magick_threads(edge_image,edge_image,edge_image->rows,1) 317 #endif 318 for (y=0; y < (ssize_t) edge_image->rows; y++) 319 { 320 register const Quantum 321 *magick_restrict p; 322 323 register ssize_t 324 x; 325 326 if (status == MagickFalse) 327 continue; 328 p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2, 329 exception); 330 if (p == (const Quantum *) NULL) 331 { 332 status=MagickFalse; 333 continue; 334 } 335 for (x=0; x < (ssize_t) edge_image->columns; x++) 336 { 337 CannyInfo 338 pixel; 339 340 double 341 dx, 342 dy; 343 344 register const Quantum 345 *magick_restrict kernel_pixels; 346 347 ssize_t 348 v; 349 350 static double 351 Gx[2][2] = 352 { 353 { -1.0, +1.0 }, 354 { -1.0, +1.0 } 355 }, 356 Gy[2][2] = 357 { 358 { +1.0, +1.0 }, 359 { -1.0, -1.0 } 360 }; 361 362 (void) ResetMagickMemory(&pixel,0,sizeof(pixel)); 363 dx=0.0; 364 dy=0.0; 365 kernel_pixels=p; 366 for (v=0; v < 2; v++) 367 { 368 ssize_t 369 u; 370 371 for (u=0; u < 2; u++) 372 { 373 double 374 intensity; 375 376 intensity=GetPixelIntensity(edge_image,kernel_pixels+u); 377 dx+=0.5*Gx[v][u]*intensity; 378 dy+=0.5*Gy[v][u]*intensity; 379 } 380 kernel_pixels+=edge_image->columns+1; 381 } 382 pixel.magnitude=hypot(dx,dy); 383 pixel.orientation=0; 384 if (fabs(dx) > MagickEpsilon) 385 { 386 double 387 slope; 388 389 slope=dy/dx; 390 if (slope < 0.0) 391 { 392 if (slope < -2.41421356237) 393 pixel.orientation=0; 394 else 395 if (slope < -0.414213562373) 396 pixel.orientation=1; 397 else 398 pixel.orientation=2; 399 } 400 else 401 { 402 if (slope > 2.41421356237) 403 pixel.orientation=0; 404 else 405 if (slope > 0.414213562373) 406 pixel.orientation=3; 407 else 408 pixel.orientation=2; 409 } 410 } 411 if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse) 412 continue; 413 p+=GetPixelChannels(edge_image); 414 } 415 } 416 edge_view=DestroyCacheView(edge_view); 417 /* 418 Non-maxima suppression, remove pixels that are not considered to be part 419 of an edge. 420 */ 421 progress=0; 422 (void) GetMatrixElement(canny_cache,0,0,&element); 423 max=element.intensity; 424 min=element.intensity; 425 edge_view=AcquireAuthenticCacheView(edge_image,exception); 426 #if defined(MAGICKCORE_OPENMP_SUPPORT) 427 #pragma omp parallel for schedule(static,4) shared(status) \ 428 magick_threads(edge_image,edge_image,edge_image->rows,1) 429 #endif 430 for (y=0; y < (ssize_t) edge_image->rows; y++) 431 { 432 register Quantum 433 *magick_restrict q; 434 435 register ssize_t 436 x; 437 438 if (status == MagickFalse) 439 continue; 440 q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1, 441 exception); 442 if (q == (Quantum *) NULL) 443 { 444 status=MagickFalse; 445 continue; 446 } 447 for (x=0; x < (ssize_t) edge_image->columns; x++) 448 { 449 CannyInfo 450 alpha_pixel, 451 beta_pixel, 452 pixel; 453 454 (void) GetMatrixElement(canny_cache,x,y,&pixel); 455 switch (pixel.orientation) 456 { 457 case 0: 458 default: 459 { 460 /* 461 0 degrees, north and south. 462 */ 463 (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel); 464 (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel); 465 break; 466 } 467 case 1: 468 { 469 /* 470 45 degrees, northwest and southeast. 471 */ 472 (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel); 473 (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel); 474 break; 475 } 476 case 2: 477 { 478 /* 479 90 degrees, east and west. 480 */ 481 (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel); 482 (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel); 483 break; 484 } 485 case 3: 486 { 487 /* 488 135 degrees, northeast and southwest. 489 */ 490 (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel); 491 (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel); 492 break; 493 } 494 } 495 pixel.intensity=pixel.magnitude; 496 if ((pixel.magnitude < alpha_pixel.magnitude) || 497 (pixel.magnitude < beta_pixel.magnitude)) 498 pixel.intensity=0; 499 (void) SetMatrixElement(canny_cache,x,y,&pixel); 500 #if defined(MAGICKCORE_OPENMP_SUPPORT) 501 #pragma omp critical (MagickCore_CannyEdgeImage) 502 #endif 503 { 504 if (pixel.intensity < min) 505 min=pixel.intensity; 506 if (pixel.intensity > max) 507 max=pixel.intensity; 508 } 509 *q=0; 510 q+=GetPixelChannels(edge_image); 511 } 512 if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse) 513 status=MagickFalse; 514 } 515 edge_view=DestroyCacheView(edge_view); 516 /* 517 Estimate hysteresis threshold. 518 */ 519 lower_threshold=lower_percent*(max-min)+min; 520 upper_threshold=upper_percent*(max-min)+min; 521 /* 522 Hysteresis threshold. 523 */ 524 edge_view=AcquireAuthenticCacheView(edge_image,exception); 525 for (y=0; y < (ssize_t) edge_image->rows; y++) 526 { 527 register ssize_t 528 x; 529 530 if (status == MagickFalse) 531 continue; 532 for (x=0; x < (ssize_t) edge_image->columns; x++) 533 { 534 CannyInfo 535 pixel; 536 537 register const Quantum 538 *magick_restrict p; 539 540 /* 541 Edge if pixel gradient higher than upper threshold. 542 */ 543 p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception); 544 if (p == (const Quantum *) NULL) 545 continue; 546 status=GetMatrixElement(canny_cache,x,y,&pixel); 547 if (status == MagickFalse) 548 continue; 549 if ((GetPixelIntensity(edge_image,p) == 0.0) && 550 (pixel.intensity >= upper_threshold)) 551 status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold, 552 exception); 553 } 554 if (image->progress_monitor != (MagickProgressMonitor) NULL) 555 { 556 MagickBooleanType 557 proceed; 558 559 #if defined(MAGICKCORE_OPENMP_SUPPORT) 560 #pragma omp critical (MagickCore_CannyEdgeImage) 561 #endif 562 proceed=SetImageProgress(image,CannyEdgeImageTag,progress++, 563 image->rows); 564 if (proceed == MagickFalse) 565 status=MagickFalse; 566 } 567 } 568 edge_view=DestroyCacheView(edge_view); 569 /* 570 Free resources. 571 */ 572 canny_cache=DestroyMatrixInfo(canny_cache); 573 return(edge_image); 574 } 575 576 /* 578 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 579 % % 580 % % 581 % % 582 % G e t I m a g e F e a t u r e s % 583 % % 584 % % 585 % % 586 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 587 % 588 % GetImageFeatures() returns features for each channel in the image in 589 % each of four directions (horizontal, vertical, left and right diagonals) 590 % for the specified distance. The features include the angular second 591 % moment, contrast, correlation, sum of squares: variance, inverse difference 592 % moment, sum average, sum varience, sum entropy, entropy, difference variance,% difference entropy, information measures of correlation 1, information 593 % measures of correlation 2, and maximum correlation coefficient. You can 594 % access the red channel contrast, for example, like this: 595 % 596 % channel_features=GetImageFeatures(image,1,exception); 597 % contrast=channel_features[RedPixelChannel].contrast[0]; 598 % 599 % Use MagickRelinquishMemory() to free the features buffer. 600 % 601 % The format of the GetImageFeatures method is: 602 % 603 % ChannelFeatures *GetImageFeatures(const Image *image, 604 % const size_t distance,ExceptionInfo *exception) 605 % 606 % A description of each parameter follows: 607 % 608 % o image: the image. 609 % 610 % o distance: the distance. 611 % 612 % o exception: return any errors or warnings in this structure. 613 % 614 */ 615 616 static inline double MagickLog10(const double x) 617 { 618 #define Log10Epsilon (1.0e-11) 619 620 if (fabs(x) < Log10Epsilon) 621 return(log10(Log10Epsilon)); 622 return(log10(fabs(x))); 623 } 624 625 MagickExport ChannelFeatures *GetImageFeatures(const Image *image, 626 const size_t distance,ExceptionInfo *exception) 627 { 628 typedef struct _ChannelStatistics 629 { 630 PixelInfo 631 direction[4]; /* horizontal, vertical, left and right diagonals */ 632 } ChannelStatistics; 633 634 CacheView 635 *image_view; 636 637 ChannelFeatures 638 *channel_features; 639 640 ChannelStatistics 641 **cooccurrence, 642 correlation, 643 *density_x, 644 *density_xy, 645 *density_y, 646 entropy_x, 647 entropy_xy, 648 entropy_xy1, 649 entropy_xy2, 650 entropy_y, 651 mean, 652 **Q, 653 *sum, 654 sum_squares, 655 variance; 656 657 PixelPacket 658 gray, 659 *grays; 660 661 MagickBooleanType 662 status; 663 664 register ssize_t 665 i, 666 r; 667 668 size_t 669 length; 670 671 unsigned int 672 number_grays; 673 674 assert(image != (Image *) NULL); 675 assert(image->signature == MagickCoreSignature); 676 if (image->debug != MagickFalse) 677 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 678 if ((image->columns < (distance+1)) || (image->rows < (distance+1))) 679 return((ChannelFeatures *) NULL); 680 length=MaxPixelChannels+1UL; 681 channel_features=(ChannelFeatures *) AcquireQuantumMemory(length, 682 sizeof(*channel_features)); 683 if (channel_features == (ChannelFeatures *) NULL) 684 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed"); 685 (void) ResetMagickMemory(channel_features,0,length* 686 sizeof(*channel_features)); 687 /* 688 Form grays. 689 */ 690 grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays)); 691 if (grays == (PixelPacket *) NULL) 692 { 693 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 694 channel_features); 695 (void) ThrowMagickException(exception,GetMagickModule(), 696 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 697 return(channel_features); 698 } 699 for (i=0; i <= (ssize_t) MaxMap; i++) 700 { 701 grays[i].red=(~0U); 702 grays[i].green=(~0U); 703 grays[i].blue=(~0U); 704 grays[i].alpha=(~0U); 705 grays[i].black=(~0U); 706 } 707 status=MagickTrue; 708 image_view=AcquireVirtualCacheView(image,exception); 709 #if defined(MAGICKCORE_OPENMP_SUPPORT) 710 #pragma omp parallel for schedule(static,4) shared(status) \ 711 magick_threads(image,image,image->rows,1) 712 #endif 713 for (r=0; r < (ssize_t) image->rows; r++) 714 { 715 register const Quantum 716 *magick_restrict p; 717 718 register ssize_t 719 x; 720 721 if (status == MagickFalse) 722 continue; 723 p=GetCacheViewVirtualPixels(image_view,0,r,image->columns,1,exception); 724 if (p == (const Quantum *) NULL) 725 { 726 status=MagickFalse; 727 continue; 728 } 729 for (x=0; x < (ssize_t) image->columns; x++) 730 { 731 grays[ScaleQuantumToMap(GetPixelRed(image,p))].red= 732 ScaleQuantumToMap(GetPixelRed(image,p)); 733 grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green= 734 ScaleQuantumToMap(GetPixelGreen(image,p)); 735 grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue= 736 ScaleQuantumToMap(GetPixelBlue(image,p)); 737 if (image->colorspace == CMYKColorspace) 738 grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black= 739 ScaleQuantumToMap(GetPixelBlack(image,p)); 740 if (image->alpha_trait != UndefinedPixelTrait) 741 grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha= 742 ScaleQuantumToMap(GetPixelAlpha(image,p)); 743 p+=GetPixelChannels(image); 744 } 745 } 746 image_view=DestroyCacheView(image_view); 747 if (status == MagickFalse) 748 { 749 grays=(PixelPacket *) RelinquishMagickMemory(grays); 750 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 751 channel_features); 752 return(channel_features); 753 } 754 (void) ResetMagickMemory(&gray,0,sizeof(gray)); 755 for (i=0; i <= (ssize_t) MaxMap; i++) 756 { 757 if (grays[i].red != ~0U) 758 grays[gray.red++].red=grays[i].red; 759 if (grays[i].green != ~0U) 760 grays[gray.green++].green=grays[i].green; 761 if (grays[i].blue != ~0U) 762 grays[gray.blue++].blue=grays[i].blue; 763 if (image->colorspace == CMYKColorspace) 764 if (grays[i].black != ~0U) 765 grays[gray.black++].black=grays[i].black; 766 if (image->alpha_trait != UndefinedPixelTrait) 767 if (grays[i].alpha != ~0U) 768 grays[gray.alpha++].alpha=grays[i].alpha; 769 } 770 /* 771 Allocate spatial dependence matrix. 772 */ 773 number_grays=gray.red; 774 if (gray.green > number_grays) 775 number_grays=gray.green; 776 if (gray.blue > number_grays) 777 number_grays=gray.blue; 778 if (image->colorspace == CMYKColorspace) 779 if (gray.black > number_grays) 780 number_grays=gray.black; 781 if (image->alpha_trait != UndefinedPixelTrait) 782 if (gray.alpha > number_grays) 783 number_grays=gray.alpha; 784 cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays, 785 sizeof(*cooccurrence)); 786 density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 787 sizeof(*density_x)); 788 density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 789 sizeof(*density_xy)); 790 density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 791 sizeof(*density_y)); 792 Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q)); 793 sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum)); 794 if ((cooccurrence == (ChannelStatistics **) NULL) || 795 (density_x == (ChannelStatistics *) NULL) || 796 (density_xy == (ChannelStatistics *) NULL) || 797 (density_y == (ChannelStatistics *) NULL) || 798 (Q == (ChannelStatistics **) NULL) || 799 (sum == (ChannelStatistics *) NULL)) 800 { 801 if (Q != (ChannelStatistics **) NULL) 802 { 803 for (i=0; i < (ssize_t) number_grays; i++) 804 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 805 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 806 } 807 if (sum != (ChannelStatistics *) NULL) 808 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 809 if (density_y != (ChannelStatistics *) NULL) 810 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 811 if (density_xy != (ChannelStatistics *) NULL) 812 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 813 if (density_x != (ChannelStatistics *) NULL) 814 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 815 if (cooccurrence != (ChannelStatistics **) NULL) 816 { 817 for (i=0; i < (ssize_t) number_grays; i++) 818 cooccurrence[i]=(ChannelStatistics *) 819 RelinquishMagickMemory(cooccurrence[i]); 820 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory( 821 cooccurrence); 822 } 823 grays=(PixelPacket *) RelinquishMagickMemory(grays); 824 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 825 channel_features); 826 (void) ThrowMagickException(exception,GetMagickModule(), 827 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 828 return(channel_features); 829 } 830 (void) ResetMagickMemory(&correlation,0,sizeof(correlation)); 831 (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x)); 832 (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy)); 833 (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y)); 834 (void) ResetMagickMemory(&mean,0,sizeof(mean)); 835 (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum)); 836 (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares)); 837 (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy)); 838 (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x)); 839 (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy)); 840 (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1)); 841 (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2)); 842 (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y)); 843 (void) ResetMagickMemory(&variance,0,sizeof(variance)); 844 for (i=0; i < (ssize_t) number_grays; i++) 845 { 846 cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays, 847 sizeof(**cooccurrence)); 848 Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q)); 849 if ((cooccurrence[i] == (ChannelStatistics *) NULL) || 850 (Q[i] == (ChannelStatistics *) NULL)) 851 break; 852 (void) ResetMagickMemory(cooccurrence[i],0,number_grays* 853 sizeof(**cooccurrence)); 854 (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q)); 855 } 856 if (i < (ssize_t) number_grays) 857 { 858 for (i--; i >= 0; i--) 859 { 860 if (Q[i] != (ChannelStatistics *) NULL) 861 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 862 if (cooccurrence[i] != (ChannelStatistics *) NULL) 863 cooccurrence[i]=(ChannelStatistics *) 864 RelinquishMagickMemory(cooccurrence[i]); 865 } 866 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 867 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 868 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 869 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 870 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 871 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 872 grays=(PixelPacket *) RelinquishMagickMemory(grays); 873 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 874 channel_features); 875 (void) ThrowMagickException(exception,GetMagickModule(), 876 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 877 return(channel_features); 878 } 879 /* 880 Initialize spatial dependence matrix. 881 */ 882 status=MagickTrue; 883 image_view=AcquireVirtualCacheView(image,exception); 884 for (r=0; r < (ssize_t) image->rows; r++) 885 { 886 register const Quantum 887 *magick_restrict p; 888 889 register ssize_t 890 x; 891 892 ssize_t 893 offset, 894 u, 895 v; 896 897 if (status == MagickFalse) 898 continue; 899 p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,r,image->columns+ 900 2*distance,distance+2,exception); 901 if (p == (const Quantum *) NULL) 902 { 903 status=MagickFalse; 904 continue; 905 } 906 p+=distance*GetPixelChannels(image);; 907 for (x=0; x < (ssize_t) image->columns; x++) 908 { 909 for (i=0; i < 4; i++) 910 { 911 switch (i) 912 { 913 case 0: 914 default: 915 { 916 /* 917 Horizontal adjacency. 918 */ 919 offset=(ssize_t) distance; 920 break; 921 } 922 case 1: 923 { 924 /* 925 Vertical adjacency. 926 */ 927 offset=(ssize_t) (image->columns+2*distance); 928 break; 929 } 930 case 2: 931 { 932 /* 933 Right diagonal adjacency. 934 */ 935 offset=(ssize_t) ((image->columns+2*distance)-distance); 936 break; 937 } 938 case 3: 939 { 940 /* 941 Left diagonal adjacency. 942 */ 943 offset=(ssize_t) ((image->columns+2*distance)+distance); 944 break; 945 } 946 } 947 u=0; 948 v=0; 949 while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p))) 950 u++; 951 while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image)))) 952 v++; 953 cooccurrence[u][v].direction[i].red++; 954 cooccurrence[v][u].direction[i].red++; 955 u=0; 956 v=0; 957 while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p))) 958 u++; 959 while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image)))) 960 v++; 961 cooccurrence[u][v].direction[i].green++; 962 cooccurrence[v][u].direction[i].green++; 963 u=0; 964 v=0; 965 while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p))) 966 u++; 967 while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image)))) 968 v++; 969 cooccurrence[u][v].direction[i].blue++; 970 cooccurrence[v][u].direction[i].blue++; 971 if (image->colorspace == CMYKColorspace) 972 { 973 u=0; 974 v=0; 975 while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p))) 976 u++; 977 while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image)))) 978 v++; 979 cooccurrence[u][v].direction[i].black++; 980 cooccurrence[v][u].direction[i].black++; 981 } 982 if (image->alpha_trait != UndefinedPixelTrait) 983 { 984 u=0; 985 v=0; 986 while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p))) 987 u++; 988 while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image)))) 989 v++; 990 cooccurrence[u][v].direction[i].alpha++; 991 cooccurrence[v][u].direction[i].alpha++; 992 } 993 } 994 p+=GetPixelChannels(image); 995 } 996 } 997 grays=(PixelPacket *) RelinquishMagickMemory(grays); 998 image_view=DestroyCacheView(image_view); 999 if (status == MagickFalse) 1000 { 1001 for (i=0; i < (ssize_t) number_grays; i++) 1002 cooccurrence[i]=(ChannelStatistics *) 1003 RelinquishMagickMemory(cooccurrence[i]); 1004 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 1005 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 1006 channel_features); 1007 (void) ThrowMagickException(exception,GetMagickModule(), 1008 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 1009 return(channel_features); 1010 } 1011 /* 1012 Normalize spatial dependence matrix. 1013 */ 1014 for (i=0; i < 4; i++) 1015 { 1016 double 1017 normalize; 1018 1019 register ssize_t 1020 y; 1021 1022 switch (i) 1023 { 1024 case 0: 1025 default: 1026 { 1027 /* 1028 Horizontal adjacency. 1029 */ 1030 normalize=2.0*image->rows*(image->columns-distance); 1031 break; 1032 } 1033 case 1: 1034 { 1035 /* 1036 Vertical adjacency. 1037 */ 1038 normalize=2.0*(image->rows-distance)*image->columns; 1039 break; 1040 } 1041 case 2: 1042 { 1043 /* 1044 Right diagonal adjacency. 1045 */ 1046 normalize=2.0*(image->rows-distance)*(image->columns-distance); 1047 break; 1048 } 1049 case 3: 1050 { 1051 /* 1052 Left diagonal adjacency. 1053 */ 1054 normalize=2.0*(image->rows-distance)*(image->columns-distance); 1055 break; 1056 } 1057 } 1058 normalize=PerceptibleReciprocal(normalize); 1059 for (y=0; y < (ssize_t) number_grays; y++) 1060 { 1061 register ssize_t 1062 x; 1063 1064 for (x=0; x < (ssize_t) number_grays; x++) 1065 { 1066 cooccurrence[x][y].direction[i].red*=normalize; 1067 cooccurrence[x][y].direction[i].green*=normalize; 1068 cooccurrence[x][y].direction[i].blue*=normalize; 1069 if (image->colorspace == CMYKColorspace) 1070 cooccurrence[x][y].direction[i].black*=normalize; 1071 if (image->alpha_trait != UndefinedPixelTrait) 1072 cooccurrence[x][y].direction[i].alpha*=normalize; 1073 } 1074 } 1075 } 1076 /* 1077 Compute texture features. 1078 */ 1079 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1080 #pragma omp parallel for schedule(static,4) shared(status) \ 1081 magick_threads(image,image,number_grays,1) 1082 #endif 1083 for (i=0; i < 4; i++) 1084 { 1085 register ssize_t 1086 y; 1087 1088 for (y=0; y < (ssize_t) number_grays; y++) 1089 { 1090 register ssize_t 1091 x; 1092 1093 for (x=0; x < (ssize_t) number_grays; x++) 1094 { 1095 /* 1096 Angular second moment: measure of homogeneity of the image. 1097 */ 1098 channel_features[RedPixelChannel].angular_second_moment[i]+= 1099 cooccurrence[x][y].direction[i].red* 1100 cooccurrence[x][y].direction[i].red; 1101 channel_features[GreenPixelChannel].angular_second_moment[i]+= 1102 cooccurrence[x][y].direction[i].green* 1103 cooccurrence[x][y].direction[i].green; 1104 channel_features[BluePixelChannel].angular_second_moment[i]+= 1105 cooccurrence[x][y].direction[i].blue* 1106 cooccurrence[x][y].direction[i].blue; 1107 if (image->colorspace == CMYKColorspace) 1108 channel_features[BlackPixelChannel].angular_second_moment[i]+= 1109 cooccurrence[x][y].direction[i].black* 1110 cooccurrence[x][y].direction[i].black; 1111 if (image->alpha_trait != UndefinedPixelTrait) 1112 channel_features[AlphaPixelChannel].angular_second_moment[i]+= 1113 cooccurrence[x][y].direction[i].alpha* 1114 cooccurrence[x][y].direction[i].alpha; 1115 /* 1116 Correlation: measure of linear-dependencies in the image. 1117 */ 1118 sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red; 1119 sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green; 1120 sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1121 if (image->colorspace == CMYKColorspace) 1122 sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black; 1123 if (image->alpha_trait != UndefinedPixelTrait) 1124 sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha; 1125 correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red; 1126 correlation.direction[i].green+=x*y* 1127 cooccurrence[x][y].direction[i].green; 1128 correlation.direction[i].blue+=x*y* 1129 cooccurrence[x][y].direction[i].blue; 1130 if (image->colorspace == CMYKColorspace) 1131 correlation.direction[i].black+=x*y* 1132 cooccurrence[x][y].direction[i].black; 1133 if (image->alpha_trait != UndefinedPixelTrait) 1134 correlation.direction[i].alpha+=x*y* 1135 cooccurrence[x][y].direction[i].alpha; 1136 /* 1137 Inverse Difference Moment. 1138 */ 1139 channel_features[RedPixelChannel].inverse_difference_moment[i]+= 1140 cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1); 1141 channel_features[GreenPixelChannel].inverse_difference_moment[i]+= 1142 cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1); 1143 channel_features[BluePixelChannel].inverse_difference_moment[i]+= 1144 cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1); 1145 if (image->colorspace == CMYKColorspace) 1146 channel_features[BlackPixelChannel].inverse_difference_moment[i]+= 1147 cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1); 1148 if (image->alpha_trait != UndefinedPixelTrait) 1149 channel_features[AlphaPixelChannel].inverse_difference_moment[i]+= 1150 cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1); 1151 /* 1152 Sum average. 1153 */ 1154 density_xy[y+x+2].direction[i].red+= 1155 cooccurrence[x][y].direction[i].red; 1156 density_xy[y+x+2].direction[i].green+= 1157 cooccurrence[x][y].direction[i].green; 1158 density_xy[y+x+2].direction[i].blue+= 1159 cooccurrence[x][y].direction[i].blue; 1160 if (image->colorspace == CMYKColorspace) 1161 density_xy[y+x+2].direction[i].black+= 1162 cooccurrence[x][y].direction[i].black; 1163 if (image->alpha_trait != UndefinedPixelTrait) 1164 density_xy[y+x+2].direction[i].alpha+= 1165 cooccurrence[x][y].direction[i].alpha; 1166 /* 1167 Entropy. 1168 */ 1169 channel_features[RedPixelChannel].entropy[i]-= 1170 cooccurrence[x][y].direction[i].red* 1171 MagickLog10(cooccurrence[x][y].direction[i].red); 1172 channel_features[GreenPixelChannel].entropy[i]-= 1173 cooccurrence[x][y].direction[i].green* 1174 MagickLog10(cooccurrence[x][y].direction[i].green); 1175 channel_features[BluePixelChannel].entropy[i]-= 1176 cooccurrence[x][y].direction[i].blue* 1177 MagickLog10(cooccurrence[x][y].direction[i].blue); 1178 if (image->colorspace == CMYKColorspace) 1179 channel_features[BlackPixelChannel].entropy[i]-= 1180 cooccurrence[x][y].direction[i].black* 1181 MagickLog10(cooccurrence[x][y].direction[i].black); 1182 if (image->alpha_trait != UndefinedPixelTrait) 1183 channel_features[AlphaPixelChannel].entropy[i]-= 1184 cooccurrence[x][y].direction[i].alpha* 1185 MagickLog10(cooccurrence[x][y].direction[i].alpha); 1186 /* 1187 Information Measures of Correlation. 1188 */ 1189 density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red; 1190 density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green; 1191 density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1192 if (image->alpha_trait != UndefinedPixelTrait) 1193 density_x[x].direction[i].alpha+= 1194 cooccurrence[x][y].direction[i].alpha; 1195 if (image->colorspace == CMYKColorspace) 1196 density_x[x].direction[i].black+= 1197 cooccurrence[x][y].direction[i].black; 1198 density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red; 1199 density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green; 1200 density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1201 if (image->colorspace == CMYKColorspace) 1202 density_y[y].direction[i].black+= 1203 cooccurrence[x][y].direction[i].black; 1204 if (image->alpha_trait != UndefinedPixelTrait) 1205 density_y[y].direction[i].alpha+= 1206 cooccurrence[x][y].direction[i].alpha; 1207 } 1208 mean.direction[i].red+=y*sum[y].direction[i].red; 1209 sum_squares.direction[i].red+=y*y*sum[y].direction[i].red; 1210 mean.direction[i].green+=y*sum[y].direction[i].green; 1211 sum_squares.direction[i].green+=y*y*sum[y].direction[i].green; 1212 mean.direction[i].blue+=y*sum[y].direction[i].blue; 1213 sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue; 1214 if (image->colorspace == CMYKColorspace) 1215 { 1216 mean.direction[i].black+=y*sum[y].direction[i].black; 1217 sum_squares.direction[i].black+=y*y*sum[y].direction[i].black; 1218 } 1219 if (image->alpha_trait != UndefinedPixelTrait) 1220 { 1221 mean.direction[i].alpha+=y*sum[y].direction[i].alpha; 1222 sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha; 1223 } 1224 } 1225 /* 1226 Correlation: measure of linear-dependencies in the image. 1227 */ 1228 channel_features[RedPixelChannel].correlation[i]= 1229 (correlation.direction[i].red-mean.direction[i].red* 1230 mean.direction[i].red)/(sqrt(sum_squares.direction[i].red- 1231 (mean.direction[i].red*mean.direction[i].red))*sqrt( 1232 sum_squares.direction[i].red-(mean.direction[i].red* 1233 mean.direction[i].red))); 1234 channel_features[GreenPixelChannel].correlation[i]= 1235 (correlation.direction[i].green-mean.direction[i].green* 1236 mean.direction[i].green)/(sqrt(sum_squares.direction[i].green- 1237 (mean.direction[i].green*mean.direction[i].green))*sqrt( 1238 sum_squares.direction[i].green-(mean.direction[i].green* 1239 mean.direction[i].green))); 1240 channel_features[BluePixelChannel].correlation[i]= 1241 (correlation.direction[i].blue-mean.direction[i].blue* 1242 mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue- 1243 (mean.direction[i].blue*mean.direction[i].blue))*sqrt( 1244 sum_squares.direction[i].blue-(mean.direction[i].blue* 1245 mean.direction[i].blue))); 1246 if (image->colorspace == CMYKColorspace) 1247 channel_features[BlackPixelChannel].correlation[i]= 1248 (correlation.direction[i].black-mean.direction[i].black* 1249 mean.direction[i].black)/(sqrt(sum_squares.direction[i].black- 1250 (mean.direction[i].black*mean.direction[i].black))*sqrt( 1251 sum_squares.direction[i].black-(mean.direction[i].black* 1252 mean.direction[i].black))); 1253 if (image->alpha_trait != UndefinedPixelTrait) 1254 channel_features[AlphaPixelChannel].correlation[i]= 1255 (correlation.direction[i].alpha-mean.direction[i].alpha* 1256 mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha- 1257 (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt( 1258 sum_squares.direction[i].alpha-(mean.direction[i].alpha* 1259 mean.direction[i].alpha))); 1260 } 1261 /* 1262 Compute more texture features. 1263 */ 1264 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1265 #pragma omp parallel for schedule(static,4) shared(status) \ 1266 magick_threads(image,image,number_grays,1) 1267 #endif 1268 for (i=0; i < 4; i++) 1269 { 1270 register ssize_t 1271 x; 1272 1273 for (x=2; x < (ssize_t) (2*number_grays); x++) 1274 { 1275 /* 1276 Sum average. 1277 */ 1278 channel_features[RedPixelChannel].sum_average[i]+= 1279 x*density_xy[x].direction[i].red; 1280 channel_features[GreenPixelChannel].sum_average[i]+= 1281 x*density_xy[x].direction[i].green; 1282 channel_features[BluePixelChannel].sum_average[i]+= 1283 x*density_xy[x].direction[i].blue; 1284 if (image->colorspace == CMYKColorspace) 1285 channel_features[BlackPixelChannel].sum_average[i]+= 1286 x*density_xy[x].direction[i].black; 1287 if (image->alpha_trait != UndefinedPixelTrait) 1288 channel_features[AlphaPixelChannel].sum_average[i]+= 1289 x*density_xy[x].direction[i].alpha; 1290 /* 1291 Sum entropy. 1292 */ 1293 channel_features[RedPixelChannel].sum_entropy[i]-= 1294 density_xy[x].direction[i].red* 1295 MagickLog10(density_xy[x].direction[i].red); 1296 channel_features[GreenPixelChannel].sum_entropy[i]-= 1297 density_xy[x].direction[i].green* 1298 MagickLog10(density_xy[x].direction[i].green); 1299 channel_features[BluePixelChannel].sum_entropy[i]-= 1300 density_xy[x].direction[i].blue* 1301 MagickLog10(density_xy[x].direction[i].blue); 1302 if (image->colorspace == CMYKColorspace) 1303 channel_features[BlackPixelChannel].sum_entropy[i]-= 1304 density_xy[x].direction[i].black* 1305 MagickLog10(density_xy[x].direction[i].black); 1306 if (image->alpha_trait != UndefinedPixelTrait) 1307 channel_features[AlphaPixelChannel].sum_entropy[i]-= 1308 density_xy[x].direction[i].alpha* 1309 MagickLog10(density_xy[x].direction[i].alpha); 1310 /* 1311 Sum variance. 1312 */ 1313 channel_features[RedPixelChannel].sum_variance[i]+= 1314 (x-channel_features[RedPixelChannel].sum_entropy[i])* 1315 (x-channel_features[RedPixelChannel].sum_entropy[i])* 1316 density_xy[x].direction[i].red; 1317 channel_features[GreenPixelChannel].sum_variance[i]+= 1318 (x-channel_features[GreenPixelChannel].sum_entropy[i])* 1319 (x-channel_features[GreenPixelChannel].sum_entropy[i])* 1320 density_xy[x].direction[i].green; 1321 channel_features[BluePixelChannel].sum_variance[i]+= 1322 (x-channel_features[BluePixelChannel].sum_entropy[i])* 1323 (x-channel_features[BluePixelChannel].sum_entropy[i])* 1324 density_xy[x].direction[i].blue; 1325 if (image->colorspace == CMYKColorspace) 1326 channel_features[BlackPixelChannel].sum_variance[i]+= 1327 (x-channel_features[BlackPixelChannel].sum_entropy[i])* 1328 (x-channel_features[BlackPixelChannel].sum_entropy[i])* 1329 density_xy[x].direction[i].black; 1330 if (image->alpha_trait != UndefinedPixelTrait) 1331 channel_features[AlphaPixelChannel].sum_variance[i]+= 1332 (x-channel_features[AlphaPixelChannel].sum_entropy[i])* 1333 (x-channel_features[AlphaPixelChannel].sum_entropy[i])* 1334 density_xy[x].direction[i].alpha; 1335 } 1336 } 1337 /* 1338 Compute more texture features. 1339 */ 1340 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1341 #pragma omp parallel for schedule(static,4) shared(status) \ 1342 magick_threads(image,image,number_grays,1) 1343 #endif 1344 for (i=0; i < 4; i++) 1345 { 1346 register ssize_t 1347 y; 1348 1349 for (y=0; y < (ssize_t) number_grays; y++) 1350 { 1351 register ssize_t 1352 x; 1353 1354 for (x=0; x < (ssize_t) number_grays; x++) 1355 { 1356 /* 1357 Sum of Squares: Variance 1358 */ 1359 variance.direction[i].red+=(y-mean.direction[i].red+1)* 1360 (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red; 1361 variance.direction[i].green+=(y-mean.direction[i].green+1)* 1362 (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green; 1363 variance.direction[i].blue+=(y-mean.direction[i].blue+1)* 1364 (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue; 1365 if (image->colorspace == CMYKColorspace) 1366 variance.direction[i].black+=(y-mean.direction[i].black+1)* 1367 (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black; 1368 if (image->alpha_trait != UndefinedPixelTrait) 1369 variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)* 1370 (y-mean.direction[i].alpha+1)* 1371 cooccurrence[x][y].direction[i].alpha; 1372 /* 1373 Sum average / Difference Variance. 1374 */ 1375 density_xy[MagickAbsoluteValue(y-x)].direction[i].red+= 1376 cooccurrence[x][y].direction[i].red; 1377 density_xy[MagickAbsoluteValue(y-x)].direction[i].green+= 1378 cooccurrence[x][y].direction[i].green; 1379 density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+= 1380 cooccurrence[x][y].direction[i].blue; 1381 if (image->colorspace == CMYKColorspace) 1382 density_xy[MagickAbsoluteValue(y-x)].direction[i].black+= 1383 cooccurrence[x][y].direction[i].black; 1384 if (image->alpha_trait != UndefinedPixelTrait) 1385 density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+= 1386 cooccurrence[x][y].direction[i].alpha; 1387 /* 1388 Information Measures of Correlation. 1389 */ 1390 entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red* 1391 MagickLog10(cooccurrence[x][y].direction[i].red); 1392 entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green* 1393 MagickLog10(cooccurrence[x][y].direction[i].green); 1394 entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue* 1395 MagickLog10(cooccurrence[x][y].direction[i].blue); 1396 if (image->colorspace == CMYKColorspace) 1397 entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black* 1398 MagickLog10(cooccurrence[x][y].direction[i].black); 1399 if (image->alpha_trait != UndefinedPixelTrait) 1400 entropy_xy.direction[i].alpha-= 1401 cooccurrence[x][y].direction[i].alpha*MagickLog10( 1402 cooccurrence[x][y].direction[i].alpha); 1403 entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red* 1404 MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red)); 1405 entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green* 1406 MagickLog10(density_x[x].direction[i].green* 1407 density_y[y].direction[i].green)); 1408 entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue* 1409 MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue)); 1410 if (image->colorspace == CMYKColorspace) 1411 entropy_xy1.direction[i].black-=( 1412 cooccurrence[x][y].direction[i].black*MagickLog10( 1413 density_x[x].direction[i].black*density_y[y].direction[i].black)); 1414 if (image->alpha_trait != UndefinedPixelTrait) 1415 entropy_xy1.direction[i].alpha-=( 1416 cooccurrence[x][y].direction[i].alpha*MagickLog10( 1417 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha)); 1418 entropy_xy2.direction[i].red-=(density_x[x].direction[i].red* 1419 density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red* 1420 density_y[y].direction[i].red)); 1421 entropy_xy2.direction[i].green-=(density_x[x].direction[i].green* 1422 density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green* 1423 density_y[y].direction[i].green)); 1424 entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue* 1425 density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue* 1426 density_y[y].direction[i].blue)); 1427 if (image->colorspace == CMYKColorspace) 1428 entropy_xy2.direction[i].black-=(density_x[x].direction[i].black* 1429 density_y[y].direction[i].black*MagickLog10( 1430 density_x[x].direction[i].black*density_y[y].direction[i].black)); 1431 if (image->alpha_trait != UndefinedPixelTrait) 1432 entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha* 1433 density_y[y].direction[i].alpha*MagickLog10( 1434 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha)); 1435 } 1436 } 1437 channel_features[RedPixelChannel].variance_sum_of_squares[i]= 1438 variance.direction[i].red; 1439 channel_features[GreenPixelChannel].variance_sum_of_squares[i]= 1440 variance.direction[i].green; 1441 channel_features[BluePixelChannel].variance_sum_of_squares[i]= 1442 variance.direction[i].blue; 1443 if (image->colorspace == CMYKColorspace) 1444 channel_features[BlackPixelChannel].variance_sum_of_squares[i]= 1445 variance.direction[i].black; 1446 if (image->alpha_trait != UndefinedPixelTrait) 1447 channel_features[AlphaPixelChannel].variance_sum_of_squares[i]= 1448 variance.direction[i].alpha; 1449 } 1450 /* 1451 Compute more texture features. 1452 */ 1453 (void) ResetMagickMemory(&variance,0,sizeof(variance)); 1454 (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares)); 1455 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1456 #pragma omp parallel for schedule(static,4) shared(status) \ 1457 magick_threads(image,image,number_grays,1) 1458 #endif 1459 for (i=0; i < 4; i++) 1460 { 1461 register ssize_t 1462 x; 1463 1464 for (x=0; x < (ssize_t) number_grays; x++) 1465 { 1466 /* 1467 Difference variance. 1468 */ 1469 variance.direction[i].red+=density_xy[x].direction[i].red; 1470 variance.direction[i].green+=density_xy[x].direction[i].green; 1471 variance.direction[i].blue+=density_xy[x].direction[i].blue; 1472 if (image->colorspace == CMYKColorspace) 1473 variance.direction[i].black+=density_xy[x].direction[i].black; 1474 if (image->alpha_trait != UndefinedPixelTrait) 1475 variance.direction[i].alpha+=density_xy[x].direction[i].alpha; 1476 sum_squares.direction[i].red+=density_xy[x].direction[i].red* 1477 density_xy[x].direction[i].red; 1478 sum_squares.direction[i].green+=density_xy[x].direction[i].green* 1479 density_xy[x].direction[i].green; 1480 sum_squares.direction[i].blue+=density_xy[x].direction[i].blue* 1481 density_xy[x].direction[i].blue; 1482 if (image->colorspace == CMYKColorspace) 1483 sum_squares.direction[i].black+=density_xy[x].direction[i].black* 1484 density_xy[x].direction[i].black; 1485 if (image->alpha_trait != UndefinedPixelTrait) 1486 sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha* 1487 density_xy[x].direction[i].alpha; 1488 /* 1489 Difference entropy. 1490 */ 1491 channel_features[RedPixelChannel].difference_entropy[i]-= 1492 density_xy[x].direction[i].red* 1493 MagickLog10(density_xy[x].direction[i].red); 1494 channel_features[GreenPixelChannel].difference_entropy[i]-= 1495 density_xy[x].direction[i].green* 1496 MagickLog10(density_xy[x].direction[i].green); 1497 channel_features[BluePixelChannel].difference_entropy[i]-= 1498 density_xy[x].direction[i].blue* 1499 MagickLog10(density_xy[x].direction[i].blue); 1500 if (image->colorspace == CMYKColorspace) 1501 channel_features[BlackPixelChannel].difference_entropy[i]-= 1502 density_xy[x].direction[i].black* 1503 MagickLog10(density_xy[x].direction[i].black); 1504 if (image->alpha_trait != UndefinedPixelTrait) 1505 channel_features[AlphaPixelChannel].difference_entropy[i]-= 1506 density_xy[x].direction[i].alpha* 1507 MagickLog10(density_xy[x].direction[i].alpha); 1508 /* 1509 Information Measures of Correlation. 1510 */ 1511 entropy_x.direction[i].red-=(density_x[x].direction[i].red* 1512 MagickLog10(density_x[x].direction[i].red)); 1513 entropy_x.direction[i].green-=(density_x[x].direction[i].green* 1514 MagickLog10(density_x[x].direction[i].green)); 1515 entropy_x.direction[i].blue-=(density_x[x].direction[i].blue* 1516 MagickLog10(density_x[x].direction[i].blue)); 1517 if (image->colorspace == CMYKColorspace) 1518 entropy_x.direction[i].black-=(density_x[x].direction[i].black* 1519 MagickLog10(density_x[x].direction[i].black)); 1520 if (image->alpha_trait != UndefinedPixelTrait) 1521 entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha* 1522 MagickLog10(density_x[x].direction[i].alpha)); 1523 entropy_y.direction[i].red-=(density_y[x].direction[i].red* 1524 MagickLog10(density_y[x].direction[i].red)); 1525 entropy_y.direction[i].green-=(density_y[x].direction[i].green* 1526 MagickLog10(density_y[x].direction[i].green)); 1527 entropy_y.direction[i].blue-=(density_y[x].direction[i].blue* 1528 MagickLog10(density_y[x].direction[i].blue)); 1529 if (image->colorspace == CMYKColorspace) 1530 entropy_y.direction[i].black-=(density_y[x].direction[i].black* 1531 MagickLog10(density_y[x].direction[i].black)); 1532 if (image->alpha_trait != UndefinedPixelTrait) 1533 entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha* 1534 MagickLog10(density_y[x].direction[i].alpha)); 1535 } 1536 /* 1537 Difference variance. 1538 */ 1539 channel_features[RedPixelChannel].difference_variance[i]= 1540 (((double) number_grays*number_grays*sum_squares.direction[i].red)- 1541 (variance.direction[i].red*variance.direction[i].red))/ 1542 ((double) number_grays*number_grays*number_grays*number_grays); 1543 channel_features[GreenPixelChannel].difference_variance[i]= 1544 (((double) number_grays*number_grays*sum_squares.direction[i].green)- 1545 (variance.direction[i].green*variance.direction[i].green))/ 1546 ((double) number_grays*number_grays*number_grays*number_grays); 1547 channel_features[BluePixelChannel].difference_variance[i]= 1548 (((double) number_grays*number_grays*sum_squares.direction[i].blue)- 1549 (variance.direction[i].blue*variance.direction[i].blue))/ 1550 ((double) number_grays*number_grays*number_grays*number_grays); 1551 if (image->colorspace == CMYKColorspace) 1552 channel_features[BlackPixelChannel].difference_variance[i]= 1553 (((double) number_grays*number_grays*sum_squares.direction[i].black)- 1554 (variance.direction[i].black*variance.direction[i].black))/ 1555 ((double) number_grays*number_grays*number_grays*number_grays); 1556 if (image->alpha_trait != UndefinedPixelTrait) 1557 channel_features[AlphaPixelChannel].difference_variance[i]= 1558 (((double) number_grays*number_grays*sum_squares.direction[i].alpha)- 1559 (variance.direction[i].alpha*variance.direction[i].alpha))/ 1560 ((double) number_grays*number_grays*number_grays*number_grays); 1561 /* 1562 Information Measures of Correlation. 1563 */ 1564 channel_features[RedPixelChannel].measure_of_correlation_1[i]= 1565 (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/ 1566 (entropy_x.direction[i].red > entropy_y.direction[i].red ? 1567 entropy_x.direction[i].red : entropy_y.direction[i].red); 1568 channel_features[GreenPixelChannel].measure_of_correlation_1[i]= 1569 (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/ 1570 (entropy_x.direction[i].green > entropy_y.direction[i].green ? 1571 entropy_x.direction[i].green : entropy_y.direction[i].green); 1572 channel_features[BluePixelChannel].measure_of_correlation_1[i]= 1573 (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/ 1574 (entropy_x.direction[i].blue > entropy_y.direction[i].blue ? 1575 entropy_x.direction[i].blue : entropy_y.direction[i].blue); 1576 if (image->colorspace == CMYKColorspace) 1577 channel_features[BlackPixelChannel].measure_of_correlation_1[i]= 1578 (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/ 1579 (entropy_x.direction[i].black > entropy_y.direction[i].black ? 1580 entropy_x.direction[i].black : entropy_y.direction[i].black); 1581 if (image->alpha_trait != UndefinedPixelTrait) 1582 channel_features[AlphaPixelChannel].measure_of_correlation_1[i]= 1583 (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/ 1584 (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ? 1585 entropy_x.direction[i].alpha : entropy_y.direction[i].alpha); 1586 channel_features[RedPixelChannel].measure_of_correlation_2[i]= 1587 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].red- 1588 entropy_xy.direction[i].red))))); 1589 channel_features[GreenPixelChannel].measure_of_correlation_2[i]= 1590 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].green- 1591 entropy_xy.direction[i].green))))); 1592 channel_features[BluePixelChannel].measure_of_correlation_2[i]= 1593 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].blue- 1594 entropy_xy.direction[i].blue))))); 1595 if (image->colorspace == CMYKColorspace) 1596 channel_features[BlackPixelChannel].measure_of_correlation_2[i]= 1597 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].black- 1598 entropy_xy.direction[i].black))))); 1599 if (image->alpha_trait != UndefinedPixelTrait) 1600 channel_features[AlphaPixelChannel].measure_of_correlation_2[i]= 1601 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].alpha- 1602 entropy_xy.direction[i].alpha))))); 1603 } 1604 /* 1605 Compute more texture features. 1606 */ 1607 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1608 #pragma omp parallel for schedule(static,4) shared(status) \ 1609 magick_threads(image,image,number_grays,1) 1610 #endif 1611 for (i=0; i < 4; i++) 1612 { 1613 ssize_t 1614 z; 1615 1616 for (z=0; z < (ssize_t) number_grays; z++) 1617 { 1618 register ssize_t 1619 y; 1620 1621 ChannelStatistics 1622 pixel; 1623 1624 (void) ResetMagickMemory(&pixel,0,sizeof(pixel)); 1625 for (y=0; y < (ssize_t) number_grays; y++) 1626 { 1627 register ssize_t 1628 x; 1629 1630 for (x=0; x < (ssize_t) number_grays; x++) 1631 { 1632 /* 1633 Contrast: amount of local variations present in an image. 1634 */ 1635 if (((y-x) == z) || ((x-y) == z)) 1636 { 1637 pixel.direction[i].red+=cooccurrence[x][y].direction[i].red; 1638 pixel.direction[i].green+=cooccurrence[x][y].direction[i].green; 1639 pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1640 if (image->colorspace == CMYKColorspace) 1641 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black; 1642 if (image->alpha_trait != UndefinedPixelTrait) 1643 pixel.direction[i].alpha+= 1644 cooccurrence[x][y].direction[i].alpha; 1645 } 1646 /* 1647 Maximum Correlation Coefficient. 1648 */ 1649 Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red* 1650 cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/ 1651 density_y[x].direction[i].red; 1652 Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green* 1653 cooccurrence[y][x].direction[i].green/ 1654 density_x[z].direction[i].green/density_y[x].direction[i].red; 1655 Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue* 1656 cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/ 1657 density_y[x].direction[i].blue; 1658 if (image->colorspace == CMYKColorspace) 1659 Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black* 1660 cooccurrence[y][x].direction[i].black/ 1661 density_x[z].direction[i].black/density_y[x].direction[i].black; 1662 if (image->alpha_trait != UndefinedPixelTrait) 1663 Q[z][y].direction[i].alpha+= 1664 cooccurrence[z][x].direction[i].alpha* 1665 cooccurrence[y][x].direction[i].alpha/ 1666 density_x[z].direction[i].alpha/ 1667 density_y[x].direction[i].alpha; 1668 } 1669 } 1670 channel_features[RedPixelChannel].contrast[i]+=z*z* 1671 pixel.direction[i].red; 1672 channel_features[GreenPixelChannel].contrast[i]+=z*z* 1673 pixel.direction[i].green; 1674 channel_features[BluePixelChannel].contrast[i]+=z*z* 1675 pixel.direction[i].blue; 1676 if (image->colorspace == CMYKColorspace) 1677 channel_features[BlackPixelChannel].contrast[i]+=z*z* 1678 pixel.direction[i].black; 1679 if (image->alpha_trait != UndefinedPixelTrait) 1680 channel_features[AlphaPixelChannel].contrast[i]+=z*z* 1681 pixel.direction[i].alpha; 1682 } 1683 /* 1684 Maximum Correlation Coefficient. 1685 Future: return second largest eigenvalue of Q. 1686 */ 1687 channel_features[RedPixelChannel].maximum_correlation_coefficient[i]= 1688 sqrt((double) -1.0); 1689 channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]= 1690 sqrt((double) -1.0); 1691 channel_features[BluePixelChannel].maximum_correlation_coefficient[i]= 1692 sqrt((double) -1.0); 1693 if (image->colorspace == CMYKColorspace) 1694 channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]= 1695 sqrt((double) -1.0); 1696 if (image->alpha_trait != UndefinedPixelTrait) 1697 channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]= 1698 sqrt((double) -1.0); 1699 } 1700 /* 1701 Relinquish resources. 1702 */ 1703 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 1704 for (i=0; i < (ssize_t) number_grays; i++) 1705 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 1706 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 1707 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 1708 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 1709 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 1710 for (i=0; i < (ssize_t) number_grays; i++) 1711 cooccurrence[i]=(ChannelStatistics *) 1712 RelinquishMagickMemory(cooccurrence[i]); 1713 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 1714 return(channel_features); 1715 } 1716 1717 /* 1719 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1720 % % 1721 % % 1722 % % 1723 % H o u g h L i n e I m a g e % 1724 % % 1725 % % 1726 % % 1727 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1728 % 1729 % Use HoughLineImage() in conjunction with any binary edge extracted image (we 1730 % recommand Canny) to identify lines in the image. The algorithm accumulates 1731 % counts for every white pixel for every possible orientation (for angles from 1732 % 0 to 179 in 1 degree increments) and distance from the center of the image to 1733 % the corner (in 1 px increments) and stores the counts in an accumulator matrix 1734 % of angle vs distance. The size of the accumulator is 180x(diagonal/2). Next 1735 % it searches this space for peaks in counts and converts the locations of the 1736 % peaks to slope and intercept in the normal x,y input image space. Use the 1737 % slope/intercepts to find the endpoints clipped to the bounds of the image. The 1738 % lines are then drawn. The counts are a measure of the length of the lines 1739 % 1740 % The format of the HoughLineImage method is: 1741 % 1742 % Image *HoughLineImage(const Image *image,const size_t width, 1743 % const size_t height,const size_t threshold,ExceptionInfo *exception) 1744 % 1745 % A description of each parameter follows: 1746 % 1747 % o image: the image. 1748 % 1749 % o width, height: find line pairs as local maxima in this neighborhood. 1750 % 1751 % o threshold: the line count threshold. 1752 % 1753 % o exception: return any errors or warnings in this structure. 1754 % 1755 */ 1756 1757 static inline double MagickRound(double x) 1758 { 1759 /* 1760 Round the fraction to nearest integer. 1761 */ 1762 if ((x-floor(x)) < (ceil(x)-x)) 1763 return(floor(x)); 1764 return(ceil(x)); 1765 } 1766 1767 static Image *RenderHoughLines(const ImageInfo *image_info,const size_t columns, 1768 const size_t rows,ExceptionInfo *exception) 1769 { 1770 #define BoundingBox "viewbox" 1771 1772 DrawInfo 1773 *draw_info; 1774 1775 Image 1776 *image; 1777 1778 MagickBooleanType 1779 status; 1780 1781 /* 1782 Open image. 1783 */ 1784 image=AcquireImage(image_info,exception); 1785 status=OpenBlob(image_info,image,ReadBinaryBlobMode,exception); 1786 if (status == MagickFalse) 1787 { 1788 image=DestroyImageList(image); 1789 return((Image *) NULL); 1790 } 1791 image->columns=columns; 1792 image->rows=rows; 1793 draw_info=CloneDrawInfo(image_info,(DrawInfo *) NULL); 1794 draw_info->affine.sx=image->resolution.x == 0.0 ? 1.0 : image->resolution.x/ 1795 DefaultResolution; 1796 draw_info->affine.sy=image->resolution.y == 0.0 ? 1.0 : image->resolution.y/ 1797 DefaultResolution; 1798 image->columns=(size_t) (draw_info->affine.sx*image->columns); 1799 image->rows=(size_t) (draw_info->affine.sy*image->rows); 1800 status=SetImageExtent(image,image->columns,image->rows,exception); 1801 if (status == MagickFalse) 1802 return(DestroyImageList(image)); 1803 if (SetImageBackgroundColor(image,exception) == MagickFalse) 1804 { 1805 image=DestroyImageList(image); 1806 return((Image *) NULL); 1807 } 1808 /* 1809 Render drawing. 1810 */ 1811 if (GetBlobStreamData(image) == (unsigned char *) NULL) 1812 draw_info->primitive=FileToString(image->filename,~0UL,exception); 1813 else 1814 { 1815 draw_info->primitive=(char *) AcquireMagickMemory((size_t) 1816 GetBlobSize(image)+1); 1817 if (draw_info->primitive != (char *) NULL) 1818 { 1819 (void) CopyMagickMemory(draw_info->primitive,GetBlobStreamData(image), 1820 (size_t) GetBlobSize(image)); 1821 draw_info->primitive[GetBlobSize(image)]='\0'; 1822 } 1823 } 1824 (void) DrawImage(image,draw_info,exception); 1825 draw_info=DestroyDrawInfo(draw_info); 1826 (void) CloseBlob(image); 1827 return(GetFirstImageInList(image)); 1828 } 1829 1830 MagickExport Image *HoughLineImage(const Image *image,const size_t width, 1831 const size_t height,const size_t threshold,ExceptionInfo *exception) 1832 { 1833 #define HoughLineImageTag "HoughLine/Image" 1834 1835 CacheView 1836 *image_view; 1837 1838 char 1839 message[MagickPathExtent], 1840 path[MagickPathExtent]; 1841 1842 const char 1843 *artifact; 1844 1845 double 1846 hough_height; 1847 1848 Image 1849 *lines_image = NULL; 1850 1851 ImageInfo 1852 *image_info; 1853 1854 int 1855 file; 1856 1857 MagickBooleanType 1858 status; 1859 1860 MagickOffsetType 1861 progress; 1862 1863 MatrixInfo 1864 *accumulator; 1865 1866 PointInfo 1867 center; 1868 1869 register ssize_t 1870 y; 1871 1872 size_t 1873 accumulator_height, 1874 accumulator_width, 1875 line_count; 1876 1877 /* 1878 Create the accumulator. 1879 */ 1880 assert(image != (const Image *) NULL); 1881 assert(image->signature == MagickCoreSignature); 1882 if (image->debug != MagickFalse) 1883 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1884 assert(exception != (ExceptionInfo *) NULL); 1885 assert(exception->signature == MagickCoreSignature); 1886 accumulator_width=180; 1887 hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ? 1888 image->rows : image->columns))/2.0); 1889 accumulator_height=(size_t) (2.0*hough_height); 1890 accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height, 1891 sizeof(double),exception); 1892 if (accumulator == (MatrixInfo *) NULL) 1893 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed"); 1894 if (NullMatrix(accumulator) == MagickFalse) 1895 { 1896 accumulator=DestroyMatrixInfo(accumulator); 1897 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed"); 1898 } 1899 /* 1900 Populate the accumulator. 1901 */ 1902 status=MagickTrue; 1903 progress=0; 1904 center.x=(double) image->columns/2.0; 1905 center.y=(double) image->rows/2.0; 1906 image_view=AcquireVirtualCacheView(image,exception); 1907 for (y=0; y < (ssize_t) image->rows; y++) 1908 { 1909 register const Quantum 1910 *magick_restrict p; 1911 1912 register ssize_t 1913 x; 1914 1915 if (status == MagickFalse) 1916 continue; 1917 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 1918 if (p == (Quantum *) NULL) 1919 { 1920 status=MagickFalse; 1921 continue; 1922 } 1923 for (x=0; x < (ssize_t) image->columns; x++) 1924 { 1925 if (GetPixelIntensity(image,p) > (QuantumRange/2.0)) 1926 { 1927 register ssize_t 1928 i; 1929 1930 for (i=0; i < 180; i++) 1931 { 1932 double 1933 count, 1934 radius; 1935 1936 radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+ 1937 (((double) y-center.y)*sin(DegreesToRadians((double) i))); 1938 (void) GetMatrixElement(accumulator,i,(ssize_t) 1939 MagickRound(radius+hough_height),&count); 1940 count++; 1941 (void) SetMatrixElement(accumulator,i,(ssize_t) 1942 MagickRound(radius+hough_height),&count); 1943 } 1944 } 1945 p+=GetPixelChannels(image); 1946 } 1947 if (image->progress_monitor != (MagickProgressMonitor) NULL) 1948 { 1949 MagickBooleanType 1950 proceed; 1951 1952 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1953 #pragma omp critical (MagickCore_CannyEdgeImage) 1954 #endif 1955 proceed=SetImageProgress(image,CannyEdgeImageTag,progress++, 1956 image->rows); 1957 if (proceed == MagickFalse) 1958 status=MagickFalse; 1959 } 1960 } 1961 image_view=DestroyCacheView(image_view); 1962 if (status == MagickFalse) 1963 { 1964 accumulator=DestroyMatrixInfo(accumulator); 1965 return((Image *) NULL); 1966 } 1967 /* 1968 Generate line segments from accumulator. 1969 */ 1970 file=AcquireUniqueFileResource(path); 1971 if (file == -1) 1972 { 1973 accumulator=DestroyMatrixInfo(accumulator); 1974 return((Image *) NULL); 1975 } 1976 (void) FormatLocaleString(message,MagickPathExtent, 1977 "# Hough line transform: %.20gx%.20g%+.20g\n",(double) width, 1978 (double) height,(double) threshold); 1979 if (write(file,message,strlen(message)) != (ssize_t) strlen(message)) 1980 status=MagickFalse; 1981 (void) FormatLocaleString(message,MagickPathExtent, 1982 "viewbox 0 0 %.20g %.20g\n",(double) image->columns,(double) image->rows); 1983 if (write(file,message,strlen(message)) != (ssize_t) strlen(message)) 1984 status=MagickFalse; 1985 line_count=image->columns > image->rows ? image->columns/4 : image->rows/4; 1986 if (threshold != 0) 1987 line_count=threshold; 1988 for (y=0; y < (ssize_t) accumulator_height; y++) 1989 { 1990 register ssize_t 1991 x; 1992 1993 for (x=0; x < (ssize_t) accumulator_width; x++) 1994 { 1995 double 1996 count; 1997 1998 (void) GetMatrixElement(accumulator,x,y,&count); 1999 if (count >= (double) line_count) 2000 { 2001 double 2002 maxima; 2003 2004 SegmentInfo 2005 line; 2006 2007 ssize_t 2008 v; 2009 2010 /* 2011 Is point a local maxima? 2012 */ 2013 maxima=count; 2014 for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++) 2015 { 2016 ssize_t 2017 u; 2018 2019 for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++) 2020 { 2021 if ((u != 0) || (v !=0)) 2022 { 2023 (void) GetMatrixElement(accumulator,x+u,y+v,&count); 2024 if (count > maxima) 2025 { 2026 maxima=count; 2027 break; 2028 } 2029 } 2030 } 2031 if (u < (ssize_t) (width/2)) 2032 break; 2033 } 2034 (void) GetMatrixElement(accumulator,x,y,&count); 2035 if (maxima > count) 2036 continue; 2037 if ((x >= 45) && (x <= 135)) 2038 { 2039 /* 2040 y = (r-x cos(t))/sin(t) 2041 */ 2042 line.x1=0.0; 2043 line.y1=((double) (y-(accumulator_height/2.0))-((line.x1- 2044 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/ 2045 sin(DegreesToRadians((double) x))+(image->rows/2.0); 2046 line.x2=(double) image->columns; 2047 line.y2=((double) (y-(accumulator_height/2.0))-((line.x2- 2048 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/ 2049 sin(DegreesToRadians((double) x))+(image->rows/2.0); 2050 } 2051 else 2052 { 2053 /* 2054 x = (r-y cos(t))/sin(t) 2055 */ 2056 line.y1=0.0; 2057 line.x1=((double) (y-(accumulator_height/2.0))-((line.y1- 2058 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/ 2059 cos(DegreesToRadians((double) x))+(image->columns/2.0); 2060 line.y2=(double) image->rows; 2061 line.x2=((double) (y-(accumulator_height/2.0))-((line.y2- 2062 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/ 2063 cos(DegreesToRadians((double) x))+(image->columns/2.0); 2064 } 2065 (void) FormatLocaleString(message,MagickPathExtent, 2066 "line %g,%g %g,%g # %g\n",line.x1,line.y1,line.x2,line.y2,maxima); 2067 if (write(file,message,strlen(message)) != (ssize_t) strlen(message)) 2068 status=MagickFalse; 2069 } 2070 } 2071 } 2072 (void) close(file); 2073 /* 2074 Render lines to image canvas. 2075 */ 2076 image_info=AcquireImageInfo(); 2077 image_info->background_color=image->background_color; 2078 (void) FormatLocaleString(image_info->filename,MagickPathExtent,"%s",path); 2079 artifact=GetImageArtifact(image,"background"); 2080 if (artifact != (const char *) NULL) 2081 (void) SetImageOption(image_info,"background",artifact); 2082 artifact=GetImageArtifact(image,"fill"); 2083 if (artifact != (const char *) NULL) 2084 (void) SetImageOption(image_info,"fill",artifact); 2085 artifact=GetImageArtifact(image,"stroke"); 2086 if (artifact != (const char *) NULL) 2087 (void) SetImageOption(image_info,"stroke",artifact); 2088 artifact=GetImageArtifact(image,"strokewidth"); 2089 if (artifact != (const char *) NULL) 2090 (void) SetImageOption(image_info,"strokewidth",artifact); 2091 lines_image=RenderHoughLines(image_info,image->columns,image->rows,exception); 2092 artifact=GetImageArtifact(image,"hough-lines:accumulator"); 2093 if ((lines_image != (Image *) NULL) && 2094 (IsStringTrue(artifact) != MagickFalse)) 2095 { 2096 Image 2097 *accumulator_image; 2098 2099 accumulator_image=MatrixToImage(accumulator,exception); 2100 if (accumulator_image != (Image *) NULL) 2101 AppendImageToList(&lines_image,accumulator_image); 2102 } 2103 /* 2104 Free resources. 2105 */ 2106 accumulator=DestroyMatrixInfo(accumulator); 2107 image_info=DestroyImageInfo(image_info); 2108 (void) RelinquishUniqueFileResource(path); 2109 return(GetFirstImageInList(lines_image)); 2110 } 2111 2112 /* 2114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2115 % % 2116 % % 2117 % % 2118 % M e a n S h i f t I m a g e % 2119 % % 2120 % % 2121 % % 2122 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2123 % 2124 % MeanShiftImage() delineate arbitrarily shaped clusters in the image. For 2125 % each pixel, it visits all the pixels in the neighborhood specified by 2126 % the window centered at the pixel and excludes those that are outside the 2127 % radius=(window-1)/2 surrounding the pixel. From those pixels, it finds those 2128 % that are within the specified color distance from the current mean, and 2129 % computes a new x,y centroid from those coordinates and a new mean. This new 2130 % x,y centroid is used as the center for a new window. This process iterates 2131 % until it converges and the final mean is replaces the (original window 2132 % center) pixel value. It repeats this process for the next pixel, etc., 2133 % until it processes all pixels in the image. Results are typically better with 2134 % colorspaces other than sRGB. We recommend YIQ, YUV or YCbCr. 2135 % 2136 % The format of the MeanShiftImage method is: 2137 % 2138 % Image *MeanShiftImage(const Image *image,const size_t width, 2139 % const size_t height,const double color_distance, 2140 % ExceptionInfo *exception) 2141 % 2142 % A description of each parameter follows: 2143 % 2144 % o image: the image. 2145 % 2146 % o width, height: find pixels in this neighborhood. 2147 % 2148 % o color_distance: the color distance. 2149 % 2150 % o exception: return any errors or warnings in this structure. 2151 % 2152 */ 2153 MagickExport Image *MeanShiftImage(const Image *image,const size_t width, 2154 const size_t height,const double color_distance,ExceptionInfo *exception) 2155 { 2156 #define MaxMeanShiftIterations 100 2157 #define MeanShiftImageTag "MeanShift/Image" 2158 2159 CacheView 2160 *image_view, 2161 *mean_view, 2162 *pixel_view; 2163 2164 Image 2165 *mean_image; 2166 2167 MagickBooleanType 2168 status; 2169 2170 MagickOffsetType 2171 progress; 2172 2173 ssize_t 2174 y; 2175 2176 assert(image != (const Image *) NULL); 2177 assert(image->signature == MagickCoreSignature); 2178 if (image->debug != MagickFalse) 2179 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 2180 assert(exception != (ExceptionInfo *) NULL); 2181 assert(exception->signature == MagickCoreSignature); 2182 mean_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception); 2183 if (mean_image == (Image *) NULL) 2184 return((Image *) NULL); 2185 if (SetImageStorageClass(mean_image,DirectClass,exception) == MagickFalse) 2186 { 2187 mean_image=DestroyImage(mean_image); 2188 return((Image *) NULL); 2189 } 2190 status=MagickTrue; 2191 progress=0; 2192 image_view=AcquireVirtualCacheView(image,exception); 2193 pixel_view=AcquireVirtualCacheView(image,exception); 2194 mean_view=AcquireAuthenticCacheView(mean_image,exception); 2195 #if defined(MAGICKCORE_OPENMP_SUPPORT) 2196 #pragma omp parallel for schedule(static,4) shared(status,progress) \ 2197 magick_threads(mean_image,mean_image,mean_image->rows,1) 2198 #endif 2199 for (y=0; y < (ssize_t) mean_image->rows; y++) 2200 { 2201 register const Quantum 2202 *magick_restrict p; 2203 2204 register Quantum 2205 *magick_restrict q; 2206 2207 register ssize_t 2208 x; 2209 2210 if (status == MagickFalse) 2211 continue; 2212 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 2213 q=GetCacheViewAuthenticPixels(mean_view,0,y,mean_image->columns,1, 2214 exception); 2215 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL)) 2216 { 2217 status=MagickFalse; 2218 continue; 2219 } 2220 for (x=0; x < (ssize_t) mean_image->columns; x++) 2221 { 2222 PixelInfo 2223 mean_pixel, 2224 previous_pixel; 2225 2226 PointInfo 2227 mean_location, 2228 previous_location; 2229 2230 register ssize_t 2231 i; 2232 2233 GetPixelInfo(image,&mean_pixel); 2234 GetPixelInfoPixel(image,p,&mean_pixel); 2235 mean_location.x=(double) x; 2236 mean_location.y=(double) y; 2237 for (i=0; i < MaxMeanShiftIterations; i++) 2238 { 2239 double 2240 distance, 2241 gamma; 2242 2243 PixelInfo 2244 sum_pixel; 2245 2246 PointInfo 2247 sum_location; 2248 2249 ssize_t 2250 count, 2251 v; 2252 2253 sum_location.x=0.0; 2254 sum_location.y=0.0; 2255 GetPixelInfo(image,&sum_pixel); 2256 previous_location=mean_location; 2257 previous_pixel=mean_pixel; 2258 count=0; 2259 for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++) 2260 { 2261 ssize_t 2262 u; 2263 2264 for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++) 2265 { 2266 if ((v*v+u*u) <= (ssize_t) ((width/2)*(height/2))) 2267 { 2268 PixelInfo 2269 pixel; 2270 2271 status=GetOneCacheViewVirtualPixelInfo(pixel_view,(ssize_t) 2272 MagickRound(mean_location.x+u),(ssize_t) MagickRound( 2273 mean_location.y+v),&pixel,exception); 2274 distance=(mean_pixel.red-pixel.red)*(mean_pixel.red-pixel.red)+ 2275 (mean_pixel.green-pixel.green)*(mean_pixel.green-pixel.green)+ 2276 (mean_pixel.blue-pixel.blue)*(mean_pixel.blue-pixel.blue); 2277 if (distance <= (color_distance*color_distance)) 2278 { 2279 sum_location.x+=mean_location.x+u; 2280 sum_location.y+=mean_location.y+v; 2281 sum_pixel.red+=pixel.red; 2282 sum_pixel.green+=pixel.green; 2283 sum_pixel.blue+=pixel.blue; 2284 sum_pixel.alpha+=pixel.alpha; 2285 count++; 2286 } 2287 } 2288 } 2289 } 2290 gamma=1.0/count; 2291 mean_location.x=gamma*sum_location.x; 2292 mean_location.y=gamma*sum_location.y; 2293 mean_pixel.red=gamma*sum_pixel.red; 2294 mean_pixel.green=gamma*sum_pixel.green; 2295 mean_pixel.blue=gamma*sum_pixel.blue; 2296 mean_pixel.alpha=gamma*sum_pixel.alpha; 2297 distance=(mean_location.x-previous_location.x)* 2298 (mean_location.x-previous_location.x)+ 2299 (mean_location.y-previous_location.y)* 2300 (mean_location.y-previous_location.y)+ 2301 255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)* 2302 255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)+ 2303 255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)* 2304 255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)+ 2305 255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue)* 2306 255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue); 2307 if (distance <= 3.0) 2308 break; 2309 } 2310 SetPixelRed(mean_image,ClampToQuantum(mean_pixel.red),q); 2311 SetPixelGreen(mean_image,ClampToQuantum(mean_pixel.green),q); 2312 SetPixelBlue(mean_image,ClampToQuantum(mean_pixel.blue),q); 2313 SetPixelAlpha(mean_image,ClampToQuantum(mean_pixel.alpha),q); 2314 p+=GetPixelChannels(image); 2315 q+=GetPixelChannels(mean_image); 2316 } 2317 if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse) 2318 status=MagickFalse; 2319 if (image->progress_monitor != (MagickProgressMonitor) NULL) 2320 { 2321 MagickBooleanType 2322 proceed; 2323 2324 #if defined(MAGICKCORE_OPENMP_SUPPORT) 2325 #pragma omp critical (MagickCore_MeanShiftImage) 2326 #endif 2327 proceed=SetImageProgress(image,MeanShiftImageTag,progress++, 2328 image->rows); 2329 if (proceed == MagickFalse) 2330 status=MagickFalse; 2331 } 2332 } 2333 mean_view=DestroyCacheView(mean_view); 2334 pixel_view=DestroyCacheView(pixel_view); 2335 image_view=DestroyCacheView(image_view); 2336 return(mean_image); 2337 } 2338
