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