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