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