feature.c revision a4898536470c21e4e82d2428acd8f9ae33dc4d9e
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% H o u g h L i n e s I m a g e % 560% % 561% % 562% % 563%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 564% 565% HoughLinesImage() identifies lines in the image. 566% 567% The format of the HoughLinesImage method is: 568% 569% Image *HoughLinesImage(const Image *image,const size_t width, 570% const size_t height,const size_t threshold,ExceptionInfo *exception) 571% 572% A description of each parameter follows: 573% 574% o image: the image. 575% 576% o width, height: find line pairs as local maxima in this neighborhood. 577% 578% o threshold: the line count threshold. 579% 580% o exception: return any errors or warnings in this structure. 581% 582*/ 583 584static inline double MagickRound(double x) 585{ 586 /* 587 Round the fraction to nearest integer. 588 */ 589 if ((x-floor(x)) < (ceil(x)-x)) 590 return(floor(x)); 591 return(ceil(x)); 592} 593 594MagickExport Image *HoughLinesImage(const Image *image,const size_t width, 595 const size_t height,const size_t threshold,ExceptionInfo *exception) 596{ 597 CacheView 598 *image_view; 599 600 char 601 message[MaxTextExtent], 602 path[MaxTextExtent]; 603 604 const char 605 *artifact; 606 607 double 608 hough_height; 609 610 Image 611 *lines_image = NULL; 612 613 ImageInfo 614 *image_info; 615 616 int 617 file; 618 619 MagickBooleanType 620 status; 621 622 MatrixInfo 623 *accumulator; 624 625 PointInfo 626 center; 627 628 register ssize_t 629 y; 630 631 size_t 632 accumulator_height, 633 accumulator_width, 634 line_count; 635 636 /* 637 Create the accumulator. 638 */ 639 assert(image != (const Image *) NULL); 640 assert(image->signature == MagickSignature); 641 if (image->debug != MagickFalse) 642 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 643 assert(exception != (ExceptionInfo *) NULL); 644 assert(exception->signature == MagickSignature); 645 accumulator_width=180; 646 hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ? 647 image->rows : image->columns))/2.0); 648 accumulator_height=(size_t) (2.0*hough_height); 649 accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height, 650 sizeof(double),exception); 651 if (accumulator == (MatrixInfo *) NULL) 652 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed"); 653 if (NullMatrix(accumulator) == MagickFalse) 654 { 655 accumulator=DestroyMatrixInfo(accumulator); 656 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed"); 657 } 658 /* 659 Populate the accumulator. 660 */ 661 status=MagickTrue; 662 center.x=(double) image->columns/2.0; 663 center.y=(double) image->rows/2.0; 664 image_view=AcquireVirtualCacheView(image,exception); 665 for (y=0; y < (ssize_t) image->rows; y++) 666 { 667 register const Quantum 668 *restrict p; 669 670 register ssize_t 671 x; 672 673 if (status == MagickFalse) 674 continue; 675 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 676 if (p == (Quantum *) NULL) 677 { 678 status=MagickFalse; 679 continue; 680 } 681 for (x=0; x < (ssize_t) image->columns; x++) 682 { 683 if (GetPixelIntensity(image,p) > (QuantumRange/2)) 684 { 685 register ssize_t 686 i; 687 688 for (i=0; i < 180; i++) 689 { 690 double 691 count, 692 radius; 693 694 radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+ 695 (((double) y-center.y)*sin(DegreesToRadians((double) i))); 696 (void) GetMatrixElement(accumulator,i,(ssize_t) 697 MagickRound(radius+hough_height),&count); 698 count++; 699 (void) SetMatrixElement(accumulator,i,(ssize_t) 700 MagickRound(radius+hough_height),&count); 701 } 702 } 703 p+=GetPixelChannels(image); 704 } 705 } 706 image_view=DestroyCacheView(image_view); 707 if (status == MagickFalse) 708 { 709 accumulator=DestroyMatrixInfo(accumulator); 710 return((Image *) NULL); 711 } 712 /* 713 Generate line segments from accumulator. 714 */ 715 file=AcquireUniqueFileResource(path); 716 if (file == -1) 717 { 718 accumulator=DestroyMatrixInfo(accumulator); 719 return((Image *) NULL); 720 } 721 (void) FormatLocaleString(message,MaxTextExtent,"viewbox 0 0 %.20g %.20g\n", 722 (double) image->columns,(double) image->rows); 723 (void) write(file,message,strlen(message)); 724 line_count=image->columns > image->rows ? image->columns/4 : image->rows/4; 725 if (threshold != 0) 726 line_count=threshold; 727 for (y=0; y < (ssize_t) accumulator_height; y++) 728 { 729 register ssize_t 730 x; 731 732 for (x=0; x < (ssize_t) accumulator_width; x++) 733 { 734 double 735 count; 736 737 (void) GetMatrixElement(accumulator,x,y,&count); 738 if (count >= (double) line_count) 739 { 740 double 741 maxima; 742 743 SegmentInfo 744 line; 745 746 ssize_t 747 v; 748 749 /* 750 Is point a local maxima? 751 */ 752 maxima=count; 753 for (v=((ssize_t) -(height/2)); v < ((ssize_t) (height/2)); v++) 754 { 755 ssize_t 756 u; 757 758 for (u=((ssize_t) -(width/2)); u < ((ssize_t) (width/2)); u++) 759 { 760 if ((u != 0) || (v !=0)) 761 { 762 (void) GetMatrixElement(accumulator,x+u,y+v,&count); 763 if (count > maxima) 764 { 765 maxima=count; 766 break; 767 } 768 } 769 } 770 if (u < (ssize_t) (width/2)) 771 break; 772 } 773 (void) GetMatrixElement(accumulator,x,y,&count); 774 if (maxima > count) 775 continue; 776 if ((x >= 45) && (x <= 135)) 777 { 778 /* 779 y = (r-x cos(t))/sin(t) 780 */ 781 line.x1=0.0; 782 line.y1=((double) (y-(accumulator_height/2.0))-((line.x1- 783 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/ 784 sin(DegreesToRadians((double) x))+(image->rows/2.0); 785 line.x2=(double) image->columns; 786 line.y2=((double) (y-(accumulator_height/2.0))-((line.x2- 787 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/ 788 sin(DegreesToRadians((double) x))+(image->rows/2.0); 789 } 790 else 791 { 792 /* 793 x = (r-y cos(t))/sin(t) 794 */ 795 line.y1=0.0; 796 line.x1=((double) (y-(accumulator_height/2.0))-((line.y1- 797 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/ 798 cos(DegreesToRadians((double) x))+(image->columns/2.0); 799 line.y2=(double) image->rows; 800 line.x2=((double) (y-(accumulator_height/2.0))-((line.y2- 801 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/ 802 cos(DegreesToRadians((double) x))+(image->columns/2.0); 803 } 804 (void) FormatLocaleString(message,MaxTextExtent, 805 "line %.20g,%.20g %.20g,%.20g # %.20g\n", 806 (double) ((ssize_t) line.x1),(double) ((ssize_t) line.y1), 807 (double) ((ssize_t) line.x2),(double) ((ssize_t) line.y2),maxima); 808 (void) write(file,message,strlen(message)); 809 } 810 } 811 } 812 (void) close(file); 813 /* 814 Render lines to image canvas. 815 */ 816 image_info=AcquireImageInfo(); 817 image_info->background_color=image->background_color; 818 (void) FormatLocaleString(image_info->filename,MaxTextExtent,"mvg:%s",path); 819 artifact=GetImageArtifact(image,"background"); 820 if (artifact != (const char *) NULL) 821 (void) SetImageOption(image_info,"background",artifact); 822 artifact=GetImageArtifact(image,"fill"); 823 if (artifact != (const char *) NULL) 824 (void) SetImageOption(image_info,"fill",artifact); 825 artifact=GetImageArtifact(image,"stroke"); 826 if (artifact != (const char *) NULL) 827 (void) SetImageOption(image_info,"stroke",artifact); 828 artifact=GetImageArtifact(image,"strokewidth"); 829 if (artifact != (const char *) NULL) 830 (void) SetImageOption(image_info,"strokewidth",artifact); 831 lines_image=ReadImage(image_info,exception); 832 artifact=GetImageArtifact(image,"hough-lines:accumulator"); 833 if ((lines_image != (Image *) NULL) && 834 (IsStringTrue(artifact) != MagickFalse)) 835 { 836 Image 837 *accumulator_image; 838 839 accumulator_image=MatrixToImage(accumulator,exception); 840 if (accumulator_image != (Image *) NULL) 841 AppendImageToList(&lines_image,accumulator_image); 842 } 843 /* 844 Free resources. 845 */ 846 accumulator=DestroyMatrixInfo(accumulator); 847 image_info=DestroyImageInfo(image_info); 848 (void) RelinquishUniqueFileResource(path); 849 return(GetFirstImageInList(lines_image)); 850} 851 852/* 853%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 854% % 855% % 856% % 857% G e t I m a g e F e a t u r e s % 858% % 859% % 860% % 861%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 862% 863% GetImageFeatures() returns features for each channel in the image in 864% each of four directions (horizontal, vertical, left and right diagonals) 865% for the specified distance. The features include the angular second 866% moment, contrast, correlation, sum of squares: variance, inverse difference 867% moment, sum average, sum varience, sum entropy, entropy, difference variance,% difference entropy, information measures of correlation 1, information 868% measures of correlation 2, and maximum correlation coefficient. You can 869% access the red channel contrast, for example, like this: 870% 871% channel_features=GetImageFeatures(image,1,exception); 872% contrast=channel_features[RedPixelChannel].contrast[0]; 873% 874% Use MagickRelinquishMemory() to free the features buffer. 875% 876% The format of the GetImageFeatures method is: 877% 878% ChannelFeatures *GetImageFeatures(const Image *image, 879% const size_t distance,ExceptionInfo *exception) 880% 881% A description of each parameter follows: 882% 883% o image: the image. 884% 885% o distance: the distance. 886% 887% o exception: return any errors or warnings in this structure. 888% 889*/ 890 891static inline ssize_t MagickAbsoluteValue(const ssize_t x) 892{ 893 if (x < 0) 894 return(-x); 895 return(x); 896} 897 898static inline double MagickLog10(const double x) 899{ 900#define Log10Epsilon (1.0e-11) 901 902 if (fabs(x) < Log10Epsilon) 903 return(log10(Log10Epsilon)); 904 return(log10(fabs(x))); 905} 906 907MagickExport ChannelFeatures *GetImageFeatures(const Image *image, 908 const size_t distance,ExceptionInfo *exception) 909{ 910 typedef struct _ChannelStatistics 911 { 912 PixelInfo 913 direction[4]; /* horizontal, vertical, left and right diagonals */ 914 } ChannelStatistics; 915 916 CacheView 917 *image_view; 918 919 ChannelFeatures 920 *channel_features; 921 922 ChannelStatistics 923 **cooccurrence, 924 correlation, 925 *density_x, 926 *density_xy, 927 *density_y, 928 entropy_x, 929 entropy_xy, 930 entropy_xy1, 931 entropy_xy2, 932 entropy_y, 933 mean, 934 **Q, 935 *sum, 936 sum_squares, 937 variance; 938 939 PixelPacket 940 gray, 941 *grays; 942 943 MagickBooleanType 944 status; 945 946 register ssize_t 947 i; 948 949 size_t 950 length; 951 952 ssize_t 953 y; 954 955 unsigned int 956 number_grays; 957 958 assert(image != (Image *) NULL); 959 assert(image->signature == MagickSignature); 960 if (image->debug != MagickFalse) 961 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 962 if ((image->columns < (distance+1)) || (image->rows < (distance+1))) 963 return((ChannelFeatures *) NULL); 964 length=CompositeChannels+1UL; 965 channel_features=(ChannelFeatures *) AcquireQuantumMemory(length, 966 sizeof(*channel_features)); 967 if (channel_features == (ChannelFeatures *) NULL) 968 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed"); 969 (void) ResetMagickMemory(channel_features,0,length* 970 sizeof(*channel_features)); 971 /* 972 Form grays. 973 */ 974 grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays)); 975 if (grays == (PixelPacket *) NULL) 976 { 977 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 978 channel_features); 979 (void) ThrowMagickException(exception,GetMagickModule(), 980 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 981 return(channel_features); 982 } 983 for (i=0; i <= (ssize_t) MaxMap; i++) 984 { 985 grays[i].red=(~0U); 986 grays[i].green=(~0U); 987 grays[i].blue=(~0U); 988 grays[i].alpha=(~0U); 989 grays[i].black=(~0U); 990 } 991 status=MagickTrue; 992 image_view=AcquireVirtualCacheView(image,exception); 993#if defined(MAGICKCORE_OPENMP_SUPPORT) 994 #pragma omp parallel for schedule(static,4) shared(status) \ 995 magick_threads(image,image,image->rows,1) 996#endif 997 for (y=0; y < (ssize_t) image->rows; y++) 998 { 999 register const Quantum 1000 *restrict p; 1001 1002 register ssize_t 1003 x; 1004 1005 if (status == MagickFalse) 1006 continue; 1007 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 1008 if (p == (const Quantum *) NULL) 1009 { 1010 status=MagickFalse; 1011 continue; 1012 } 1013 for (x=0; x < (ssize_t) image->columns; x++) 1014 { 1015 grays[ScaleQuantumToMap(GetPixelRed(image,p))].red= 1016 ScaleQuantumToMap(GetPixelRed(image,p)); 1017 grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green= 1018 ScaleQuantumToMap(GetPixelGreen(image,p)); 1019 grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue= 1020 ScaleQuantumToMap(GetPixelBlue(image,p)); 1021 if (image->colorspace == CMYKColorspace) 1022 grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black= 1023 ScaleQuantumToMap(GetPixelBlack(image,p)); 1024 if (image->alpha_trait == BlendPixelTrait) 1025 grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha= 1026 ScaleQuantumToMap(GetPixelAlpha(image,p)); 1027 p+=GetPixelChannels(image); 1028 } 1029 } 1030 image_view=DestroyCacheView(image_view); 1031 if (status == MagickFalse) 1032 { 1033 grays=(PixelPacket *) RelinquishMagickMemory(grays); 1034 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 1035 channel_features); 1036 return(channel_features); 1037 } 1038 (void) ResetMagickMemory(&gray,0,sizeof(gray)); 1039 for (i=0; i <= (ssize_t) MaxMap; i++) 1040 { 1041 if (grays[i].red != ~0U) 1042 grays[gray.red++].red=grays[i].red; 1043 if (grays[i].green != ~0U) 1044 grays[gray.green++].green=grays[i].green; 1045 if (grays[i].blue != ~0U) 1046 grays[gray.blue++].blue=grays[i].blue; 1047 if (image->colorspace == CMYKColorspace) 1048 if (grays[i].black != ~0U) 1049 grays[gray.black++].black=grays[i].black; 1050 if (image->alpha_trait == BlendPixelTrait) 1051 if (grays[i].alpha != ~0U) 1052 grays[gray.alpha++].alpha=grays[i].alpha; 1053 } 1054 /* 1055 Allocate spatial dependence matrix. 1056 */ 1057 number_grays=gray.red; 1058 if (gray.green > number_grays) 1059 number_grays=gray.green; 1060 if (gray.blue > number_grays) 1061 number_grays=gray.blue; 1062 if (image->colorspace == CMYKColorspace) 1063 if (gray.black > number_grays) 1064 number_grays=gray.black; 1065 if (image->alpha_trait == BlendPixelTrait) 1066 if (gray.alpha > number_grays) 1067 number_grays=gray.alpha; 1068 cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays, 1069 sizeof(*cooccurrence)); 1070 density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 1071 sizeof(*density_x)); 1072 density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 1073 sizeof(*density_xy)); 1074 density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1), 1075 sizeof(*density_y)); 1076 Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q)); 1077 sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum)); 1078 if ((cooccurrence == (ChannelStatistics **) NULL) || 1079 (density_x == (ChannelStatistics *) NULL) || 1080 (density_xy == (ChannelStatistics *) NULL) || 1081 (density_y == (ChannelStatistics *) NULL) || 1082 (Q == (ChannelStatistics **) NULL) || 1083 (sum == (ChannelStatistics *) NULL)) 1084 { 1085 if (Q != (ChannelStatistics **) NULL) 1086 { 1087 for (i=0; i < (ssize_t) number_grays; i++) 1088 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 1089 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 1090 } 1091 if (sum != (ChannelStatistics *) NULL) 1092 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 1093 if (density_y != (ChannelStatistics *) NULL) 1094 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 1095 if (density_xy != (ChannelStatistics *) NULL) 1096 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 1097 if (density_x != (ChannelStatistics *) NULL) 1098 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 1099 if (cooccurrence != (ChannelStatistics **) NULL) 1100 { 1101 for (i=0; i < (ssize_t) number_grays; i++) 1102 cooccurrence[i]=(ChannelStatistics *) 1103 RelinquishMagickMemory(cooccurrence[i]); 1104 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory( 1105 cooccurrence); 1106 } 1107 grays=(PixelPacket *) RelinquishMagickMemory(grays); 1108 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 1109 channel_features); 1110 (void) ThrowMagickException(exception,GetMagickModule(), 1111 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 1112 return(channel_features); 1113 } 1114 (void) ResetMagickMemory(&correlation,0,sizeof(correlation)); 1115 (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x)); 1116 (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy)); 1117 (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y)); 1118 (void) ResetMagickMemory(&mean,0,sizeof(mean)); 1119 (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum)); 1120 (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares)); 1121 (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy)); 1122 (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x)); 1123 (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy)); 1124 (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1)); 1125 (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2)); 1126 (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y)); 1127 (void) ResetMagickMemory(&variance,0,sizeof(variance)); 1128 for (i=0; i < (ssize_t) number_grays; i++) 1129 { 1130 cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays, 1131 sizeof(**cooccurrence)); 1132 Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q)); 1133 if ((cooccurrence[i] == (ChannelStatistics *) NULL) || 1134 (Q[i] == (ChannelStatistics *) NULL)) 1135 break; 1136 (void) ResetMagickMemory(cooccurrence[i],0,number_grays* 1137 sizeof(**cooccurrence)); 1138 (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q)); 1139 } 1140 if (i < (ssize_t) number_grays) 1141 { 1142 for (i--; i >= 0; i--) 1143 { 1144 if (Q[i] != (ChannelStatistics *) NULL) 1145 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 1146 if (cooccurrence[i] != (ChannelStatistics *) NULL) 1147 cooccurrence[i]=(ChannelStatistics *) 1148 RelinquishMagickMemory(cooccurrence[i]); 1149 } 1150 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 1151 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 1152 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 1153 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 1154 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 1155 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 1156 grays=(PixelPacket *) RelinquishMagickMemory(grays); 1157 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 1158 channel_features); 1159 (void) ThrowMagickException(exception,GetMagickModule(), 1160 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 1161 return(channel_features); 1162 } 1163 /* 1164 Initialize spatial dependence matrix. 1165 */ 1166 status=MagickTrue; 1167 image_view=AcquireVirtualCacheView(image,exception); 1168 for (y=0; y < (ssize_t) image->rows; y++) 1169 { 1170 register const Quantum 1171 *restrict p; 1172 1173 register ssize_t 1174 x; 1175 1176 ssize_t 1177 i, 1178 offset, 1179 u, 1180 v; 1181 1182 if (status == MagickFalse) 1183 continue; 1184 p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+ 1185 2*distance,distance+2,exception); 1186 if (p == (const Quantum *) NULL) 1187 { 1188 status=MagickFalse; 1189 continue; 1190 } 1191 p+=distance*GetPixelChannels(image);; 1192 for (x=0; x < (ssize_t) image->columns; x++) 1193 { 1194 for (i=0; i < 4; i++) 1195 { 1196 switch (i) 1197 { 1198 case 0: 1199 default: 1200 { 1201 /* 1202 Horizontal adjacency. 1203 */ 1204 offset=(ssize_t) distance; 1205 break; 1206 } 1207 case 1: 1208 { 1209 /* 1210 Vertical adjacency. 1211 */ 1212 offset=(ssize_t) (image->columns+2*distance); 1213 break; 1214 } 1215 case 2: 1216 { 1217 /* 1218 Right diagonal adjacency. 1219 */ 1220 offset=(ssize_t) ((image->columns+2*distance)-distance); 1221 break; 1222 } 1223 case 3: 1224 { 1225 /* 1226 Left diagonal adjacency. 1227 */ 1228 offset=(ssize_t) ((image->columns+2*distance)+distance); 1229 break; 1230 } 1231 } 1232 u=0; 1233 v=0; 1234 while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p))) 1235 u++; 1236 while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image)))) 1237 v++; 1238 cooccurrence[u][v].direction[i].red++; 1239 cooccurrence[v][u].direction[i].red++; 1240 u=0; 1241 v=0; 1242 while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p))) 1243 u++; 1244 while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image)))) 1245 v++; 1246 cooccurrence[u][v].direction[i].green++; 1247 cooccurrence[v][u].direction[i].green++; 1248 u=0; 1249 v=0; 1250 while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p))) 1251 u++; 1252 while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image)))) 1253 v++; 1254 cooccurrence[u][v].direction[i].blue++; 1255 cooccurrence[v][u].direction[i].blue++; 1256 if (image->colorspace == CMYKColorspace) 1257 { 1258 u=0; 1259 v=0; 1260 while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p))) 1261 u++; 1262 while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image)))) 1263 v++; 1264 cooccurrence[u][v].direction[i].black++; 1265 cooccurrence[v][u].direction[i].black++; 1266 } 1267 if (image->alpha_trait == BlendPixelTrait) 1268 { 1269 u=0; 1270 v=0; 1271 while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p))) 1272 u++; 1273 while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image)))) 1274 v++; 1275 cooccurrence[u][v].direction[i].alpha++; 1276 cooccurrence[v][u].direction[i].alpha++; 1277 } 1278 } 1279 p+=GetPixelChannels(image); 1280 } 1281 } 1282 grays=(PixelPacket *) RelinquishMagickMemory(grays); 1283 image_view=DestroyCacheView(image_view); 1284 if (status == MagickFalse) 1285 { 1286 for (i=0; i < (ssize_t) number_grays; i++) 1287 cooccurrence[i]=(ChannelStatistics *) 1288 RelinquishMagickMemory(cooccurrence[i]); 1289 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 1290 channel_features=(ChannelFeatures *) RelinquishMagickMemory( 1291 channel_features); 1292 (void) ThrowMagickException(exception,GetMagickModule(), 1293 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 1294 return(channel_features); 1295 } 1296 /* 1297 Normalize spatial dependence matrix. 1298 */ 1299 for (i=0; i < 4; i++) 1300 { 1301 double 1302 normalize; 1303 1304 register ssize_t 1305 y; 1306 1307 switch (i) 1308 { 1309 case 0: 1310 default: 1311 { 1312 /* 1313 Horizontal adjacency. 1314 */ 1315 normalize=2.0*image->rows*(image->columns-distance); 1316 break; 1317 } 1318 case 1: 1319 { 1320 /* 1321 Vertical adjacency. 1322 */ 1323 normalize=2.0*(image->rows-distance)*image->columns; 1324 break; 1325 } 1326 case 2: 1327 { 1328 /* 1329 Right diagonal adjacency. 1330 */ 1331 normalize=2.0*(image->rows-distance)*(image->columns-distance); 1332 break; 1333 } 1334 case 3: 1335 { 1336 /* 1337 Left diagonal adjacency. 1338 */ 1339 normalize=2.0*(image->rows-distance)*(image->columns-distance); 1340 break; 1341 } 1342 } 1343 normalize=PerceptibleReciprocal(normalize); 1344 for (y=0; y < (ssize_t) number_grays; y++) 1345 { 1346 register ssize_t 1347 x; 1348 1349 for (x=0; x < (ssize_t) number_grays; x++) 1350 { 1351 cooccurrence[x][y].direction[i].red*=normalize; 1352 cooccurrence[x][y].direction[i].green*=normalize; 1353 cooccurrence[x][y].direction[i].blue*=normalize; 1354 if (image->colorspace == CMYKColorspace) 1355 cooccurrence[x][y].direction[i].black*=normalize; 1356 if (image->alpha_trait == BlendPixelTrait) 1357 cooccurrence[x][y].direction[i].alpha*=normalize; 1358 } 1359 } 1360 } 1361 /* 1362 Compute texture features. 1363 */ 1364#if defined(MAGICKCORE_OPENMP_SUPPORT) 1365 #pragma omp parallel for schedule(static,4) shared(status) \ 1366 magick_threads(image,image,number_grays,1) 1367#endif 1368 for (i=0; i < 4; i++) 1369 { 1370 register ssize_t 1371 y; 1372 1373 for (y=0; y < (ssize_t) number_grays; y++) 1374 { 1375 register ssize_t 1376 x; 1377 1378 for (x=0; x < (ssize_t) number_grays; x++) 1379 { 1380 /* 1381 Angular second moment: measure of homogeneity of the image. 1382 */ 1383 channel_features[RedPixelChannel].angular_second_moment[i]+= 1384 cooccurrence[x][y].direction[i].red* 1385 cooccurrence[x][y].direction[i].red; 1386 channel_features[GreenPixelChannel].angular_second_moment[i]+= 1387 cooccurrence[x][y].direction[i].green* 1388 cooccurrence[x][y].direction[i].green; 1389 channel_features[BluePixelChannel].angular_second_moment[i]+= 1390 cooccurrence[x][y].direction[i].blue* 1391 cooccurrence[x][y].direction[i].blue; 1392 if (image->colorspace == CMYKColorspace) 1393 channel_features[BlackPixelChannel].angular_second_moment[i]+= 1394 cooccurrence[x][y].direction[i].black* 1395 cooccurrence[x][y].direction[i].black; 1396 if (image->alpha_trait == BlendPixelTrait) 1397 channel_features[AlphaPixelChannel].angular_second_moment[i]+= 1398 cooccurrence[x][y].direction[i].alpha* 1399 cooccurrence[x][y].direction[i].alpha; 1400 /* 1401 Correlation: measure of linear-dependencies in the image. 1402 */ 1403 sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red; 1404 sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green; 1405 sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1406 if (image->colorspace == CMYKColorspace) 1407 sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black; 1408 if (image->alpha_trait == BlendPixelTrait) 1409 sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha; 1410 correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red; 1411 correlation.direction[i].green+=x*y* 1412 cooccurrence[x][y].direction[i].green; 1413 correlation.direction[i].blue+=x*y* 1414 cooccurrence[x][y].direction[i].blue; 1415 if (image->colorspace == CMYKColorspace) 1416 correlation.direction[i].black+=x*y* 1417 cooccurrence[x][y].direction[i].black; 1418 if (image->alpha_trait == BlendPixelTrait) 1419 correlation.direction[i].alpha+=x*y* 1420 cooccurrence[x][y].direction[i].alpha; 1421 /* 1422 Inverse Difference Moment. 1423 */ 1424 channel_features[RedPixelChannel].inverse_difference_moment[i]+= 1425 cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1); 1426 channel_features[GreenPixelChannel].inverse_difference_moment[i]+= 1427 cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1); 1428 channel_features[BluePixelChannel].inverse_difference_moment[i]+= 1429 cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1); 1430 if (image->colorspace == CMYKColorspace) 1431 channel_features[BlackPixelChannel].inverse_difference_moment[i]+= 1432 cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1); 1433 if (image->alpha_trait == BlendPixelTrait) 1434 channel_features[AlphaPixelChannel].inverse_difference_moment[i]+= 1435 cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1); 1436 /* 1437 Sum average. 1438 */ 1439 density_xy[y+x+2].direction[i].red+= 1440 cooccurrence[x][y].direction[i].red; 1441 density_xy[y+x+2].direction[i].green+= 1442 cooccurrence[x][y].direction[i].green; 1443 density_xy[y+x+2].direction[i].blue+= 1444 cooccurrence[x][y].direction[i].blue; 1445 if (image->colorspace == CMYKColorspace) 1446 density_xy[y+x+2].direction[i].black+= 1447 cooccurrence[x][y].direction[i].black; 1448 if (image->alpha_trait == BlendPixelTrait) 1449 density_xy[y+x+2].direction[i].alpha+= 1450 cooccurrence[x][y].direction[i].alpha; 1451 /* 1452 Entropy. 1453 */ 1454 channel_features[RedPixelChannel].entropy[i]-= 1455 cooccurrence[x][y].direction[i].red* 1456 MagickLog10(cooccurrence[x][y].direction[i].red); 1457 channel_features[GreenPixelChannel].entropy[i]-= 1458 cooccurrence[x][y].direction[i].green* 1459 MagickLog10(cooccurrence[x][y].direction[i].green); 1460 channel_features[BluePixelChannel].entropy[i]-= 1461 cooccurrence[x][y].direction[i].blue* 1462 MagickLog10(cooccurrence[x][y].direction[i].blue); 1463 if (image->colorspace == CMYKColorspace) 1464 channel_features[BlackPixelChannel].entropy[i]-= 1465 cooccurrence[x][y].direction[i].black* 1466 MagickLog10(cooccurrence[x][y].direction[i].black); 1467 if (image->alpha_trait == BlendPixelTrait) 1468 channel_features[AlphaPixelChannel].entropy[i]-= 1469 cooccurrence[x][y].direction[i].alpha* 1470 MagickLog10(cooccurrence[x][y].direction[i].alpha); 1471 /* 1472 Information Measures of Correlation. 1473 */ 1474 density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red; 1475 density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green; 1476 density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1477 if (image->alpha_trait == BlendPixelTrait) 1478 density_x[x].direction[i].alpha+= 1479 cooccurrence[x][y].direction[i].alpha; 1480 if (image->colorspace == CMYKColorspace) 1481 density_x[x].direction[i].black+= 1482 cooccurrence[x][y].direction[i].black; 1483 density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red; 1484 density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green; 1485 density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1486 if (image->colorspace == CMYKColorspace) 1487 density_y[y].direction[i].black+= 1488 cooccurrence[x][y].direction[i].black; 1489 if (image->alpha_trait == BlendPixelTrait) 1490 density_y[y].direction[i].alpha+= 1491 cooccurrence[x][y].direction[i].alpha; 1492 } 1493 mean.direction[i].red+=y*sum[y].direction[i].red; 1494 sum_squares.direction[i].red+=y*y*sum[y].direction[i].red; 1495 mean.direction[i].green+=y*sum[y].direction[i].green; 1496 sum_squares.direction[i].green+=y*y*sum[y].direction[i].green; 1497 mean.direction[i].blue+=y*sum[y].direction[i].blue; 1498 sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue; 1499 if (image->colorspace == CMYKColorspace) 1500 { 1501 mean.direction[i].black+=y*sum[y].direction[i].black; 1502 sum_squares.direction[i].black+=y*y*sum[y].direction[i].black; 1503 } 1504 if (image->alpha_trait == BlendPixelTrait) 1505 { 1506 mean.direction[i].alpha+=y*sum[y].direction[i].alpha; 1507 sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha; 1508 } 1509 } 1510 /* 1511 Correlation: measure of linear-dependencies in the image. 1512 */ 1513 channel_features[RedPixelChannel].correlation[i]= 1514 (correlation.direction[i].red-mean.direction[i].red* 1515 mean.direction[i].red)/(sqrt(sum_squares.direction[i].red- 1516 (mean.direction[i].red*mean.direction[i].red))*sqrt( 1517 sum_squares.direction[i].red-(mean.direction[i].red* 1518 mean.direction[i].red))); 1519 channel_features[GreenPixelChannel].correlation[i]= 1520 (correlation.direction[i].green-mean.direction[i].green* 1521 mean.direction[i].green)/(sqrt(sum_squares.direction[i].green- 1522 (mean.direction[i].green*mean.direction[i].green))*sqrt( 1523 sum_squares.direction[i].green-(mean.direction[i].green* 1524 mean.direction[i].green))); 1525 channel_features[BluePixelChannel].correlation[i]= 1526 (correlation.direction[i].blue-mean.direction[i].blue* 1527 mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue- 1528 (mean.direction[i].blue*mean.direction[i].blue))*sqrt( 1529 sum_squares.direction[i].blue-(mean.direction[i].blue* 1530 mean.direction[i].blue))); 1531 if (image->colorspace == CMYKColorspace) 1532 channel_features[BlackPixelChannel].correlation[i]= 1533 (correlation.direction[i].black-mean.direction[i].black* 1534 mean.direction[i].black)/(sqrt(sum_squares.direction[i].black- 1535 (mean.direction[i].black*mean.direction[i].black))*sqrt( 1536 sum_squares.direction[i].black-(mean.direction[i].black* 1537 mean.direction[i].black))); 1538 if (image->alpha_trait == BlendPixelTrait) 1539 channel_features[AlphaPixelChannel].correlation[i]= 1540 (correlation.direction[i].alpha-mean.direction[i].alpha* 1541 mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha- 1542 (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt( 1543 sum_squares.direction[i].alpha-(mean.direction[i].alpha* 1544 mean.direction[i].alpha))); 1545 } 1546 /* 1547 Compute more texture features. 1548 */ 1549#if defined(MAGICKCORE_OPENMP_SUPPORT) 1550 #pragma omp parallel for schedule(static,4) shared(status) \ 1551 magick_threads(image,image,number_grays,1) 1552#endif 1553 for (i=0; i < 4; i++) 1554 { 1555 register ssize_t 1556 x; 1557 1558 for (x=2; x < (ssize_t) (2*number_grays); x++) 1559 { 1560 /* 1561 Sum average. 1562 */ 1563 channel_features[RedPixelChannel].sum_average[i]+= 1564 x*density_xy[x].direction[i].red; 1565 channel_features[GreenPixelChannel].sum_average[i]+= 1566 x*density_xy[x].direction[i].green; 1567 channel_features[BluePixelChannel].sum_average[i]+= 1568 x*density_xy[x].direction[i].blue; 1569 if (image->colorspace == CMYKColorspace) 1570 channel_features[BlackPixelChannel].sum_average[i]+= 1571 x*density_xy[x].direction[i].black; 1572 if (image->alpha_trait == BlendPixelTrait) 1573 channel_features[AlphaPixelChannel].sum_average[i]+= 1574 x*density_xy[x].direction[i].alpha; 1575 /* 1576 Sum entropy. 1577 */ 1578 channel_features[RedPixelChannel].sum_entropy[i]-= 1579 density_xy[x].direction[i].red* 1580 MagickLog10(density_xy[x].direction[i].red); 1581 channel_features[GreenPixelChannel].sum_entropy[i]-= 1582 density_xy[x].direction[i].green* 1583 MagickLog10(density_xy[x].direction[i].green); 1584 channel_features[BluePixelChannel].sum_entropy[i]-= 1585 density_xy[x].direction[i].blue* 1586 MagickLog10(density_xy[x].direction[i].blue); 1587 if (image->colorspace == CMYKColorspace) 1588 channel_features[BlackPixelChannel].sum_entropy[i]-= 1589 density_xy[x].direction[i].black* 1590 MagickLog10(density_xy[x].direction[i].black); 1591 if (image->alpha_trait == BlendPixelTrait) 1592 channel_features[AlphaPixelChannel].sum_entropy[i]-= 1593 density_xy[x].direction[i].alpha* 1594 MagickLog10(density_xy[x].direction[i].alpha); 1595 /* 1596 Sum variance. 1597 */ 1598 channel_features[RedPixelChannel].sum_variance[i]+= 1599 (x-channel_features[RedPixelChannel].sum_entropy[i])* 1600 (x-channel_features[RedPixelChannel].sum_entropy[i])* 1601 density_xy[x].direction[i].red; 1602 channel_features[GreenPixelChannel].sum_variance[i]+= 1603 (x-channel_features[GreenPixelChannel].sum_entropy[i])* 1604 (x-channel_features[GreenPixelChannel].sum_entropy[i])* 1605 density_xy[x].direction[i].green; 1606 channel_features[BluePixelChannel].sum_variance[i]+= 1607 (x-channel_features[BluePixelChannel].sum_entropy[i])* 1608 (x-channel_features[BluePixelChannel].sum_entropy[i])* 1609 density_xy[x].direction[i].blue; 1610 if (image->colorspace == CMYKColorspace) 1611 channel_features[BlackPixelChannel].sum_variance[i]+= 1612 (x-channel_features[BlackPixelChannel].sum_entropy[i])* 1613 (x-channel_features[BlackPixelChannel].sum_entropy[i])* 1614 density_xy[x].direction[i].black; 1615 if (image->alpha_trait == BlendPixelTrait) 1616 channel_features[AlphaPixelChannel].sum_variance[i]+= 1617 (x-channel_features[AlphaPixelChannel].sum_entropy[i])* 1618 (x-channel_features[AlphaPixelChannel].sum_entropy[i])* 1619 density_xy[x].direction[i].alpha; 1620 } 1621 } 1622 /* 1623 Compute more texture features. 1624 */ 1625#if defined(MAGICKCORE_OPENMP_SUPPORT) 1626 #pragma omp parallel for schedule(static,4) shared(status) \ 1627 magick_threads(image,image,number_grays,1) 1628#endif 1629 for (i=0; i < 4; i++) 1630 { 1631 register ssize_t 1632 y; 1633 1634 for (y=0; y < (ssize_t) number_grays; y++) 1635 { 1636 register ssize_t 1637 x; 1638 1639 for (x=0; x < (ssize_t) number_grays; x++) 1640 { 1641 /* 1642 Sum of Squares: Variance 1643 */ 1644 variance.direction[i].red+=(y-mean.direction[i].red+1)* 1645 (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red; 1646 variance.direction[i].green+=(y-mean.direction[i].green+1)* 1647 (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green; 1648 variance.direction[i].blue+=(y-mean.direction[i].blue+1)* 1649 (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue; 1650 if (image->colorspace == CMYKColorspace) 1651 variance.direction[i].black+=(y-mean.direction[i].black+1)* 1652 (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black; 1653 if (image->alpha_trait == BlendPixelTrait) 1654 variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)* 1655 (y-mean.direction[i].alpha+1)* 1656 cooccurrence[x][y].direction[i].alpha; 1657 /* 1658 Sum average / Difference Variance. 1659 */ 1660 density_xy[MagickAbsoluteValue(y-x)].direction[i].red+= 1661 cooccurrence[x][y].direction[i].red; 1662 density_xy[MagickAbsoluteValue(y-x)].direction[i].green+= 1663 cooccurrence[x][y].direction[i].green; 1664 density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+= 1665 cooccurrence[x][y].direction[i].blue; 1666 if (image->colorspace == CMYKColorspace) 1667 density_xy[MagickAbsoluteValue(y-x)].direction[i].black+= 1668 cooccurrence[x][y].direction[i].black; 1669 if (image->alpha_trait == BlendPixelTrait) 1670 density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+= 1671 cooccurrence[x][y].direction[i].alpha; 1672 /* 1673 Information Measures of Correlation. 1674 */ 1675 entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red* 1676 MagickLog10(cooccurrence[x][y].direction[i].red); 1677 entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green* 1678 MagickLog10(cooccurrence[x][y].direction[i].green); 1679 entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue* 1680 MagickLog10(cooccurrence[x][y].direction[i].blue); 1681 if (image->colorspace == CMYKColorspace) 1682 entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black* 1683 MagickLog10(cooccurrence[x][y].direction[i].black); 1684 if (image->alpha_trait == BlendPixelTrait) 1685 entropy_xy.direction[i].alpha-= 1686 cooccurrence[x][y].direction[i].alpha*MagickLog10( 1687 cooccurrence[x][y].direction[i].alpha); 1688 entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red* 1689 MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red)); 1690 entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green* 1691 MagickLog10(density_x[x].direction[i].green* 1692 density_y[y].direction[i].green)); 1693 entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue* 1694 MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue)); 1695 if (image->colorspace == CMYKColorspace) 1696 entropy_xy1.direction[i].black-=( 1697 cooccurrence[x][y].direction[i].black*MagickLog10( 1698 density_x[x].direction[i].black*density_y[y].direction[i].black)); 1699 if (image->alpha_trait == BlendPixelTrait) 1700 entropy_xy1.direction[i].alpha-=( 1701 cooccurrence[x][y].direction[i].alpha*MagickLog10( 1702 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha)); 1703 entropy_xy2.direction[i].red-=(density_x[x].direction[i].red* 1704 density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red* 1705 density_y[y].direction[i].red)); 1706 entropy_xy2.direction[i].green-=(density_x[x].direction[i].green* 1707 density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green* 1708 density_y[y].direction[i].green)); 1709 entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue* 1710 density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue* 1711 density_y[y].direction[i].blue)); 1712 if (image->colorspace == CMYKColorspace) 1713 entropy_xy2.direction[i].black-=(density_x[x].direction[i].black* 1714 density_y[y].direction[i].black*MagickLog10( 1715 density_x[x].direction[i].black*density_y[y].direction[i].black)); 1716 if (image->alpha_trait == BlendPixelTrait) 1717 entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha* 1718 density_y[y].direction[i].alpha*MagickLog10( 1719 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha)); 1720 } 1721 } 1722 channel_features[RedPixelChannel].variance_sum_of_squares[i]= 1723 variance.direction[i].red; 1724 channel_features[GreenPixelChannel].variance_sum_of_squares[i]= 1725 variance.direction[i].green; 1726 channel_features[BluePixelChannel].variance_sum_of_squares[i]= 1727 variance.direction[i].blue; 1728 if (image->colorspace == CMYKColorspace) 1729 channel_features[BlackPixelChannel].variance_sum_of_squares[i]= 1730 variance.direction[i].black; 1731 if (image->alpha_trait == BlendPixelTrait) 1732 channel_features[AlphaPixelChannel].variance_sum_of_squares[i]= 1733 variance.direction[i].alpha; 1734 } 1735 /* 1736 Compute more texture features. 1737 */ 1738 (void) ResetMagickMemory(&variance,0,sizeof(variance)); 1739 (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares)); 1740#if defined(MAGICKCORE_OPENMP_SUPPORT) 1741 #pragma omp parallel for schedule(static,4) shared(status) \ 1742 magick_threads(image,image,number_grays,1) 1743#endif 1744 for (i=0; i < 4; i++) 1745 { 1746 register ssize_t 1747 x; 1748 1749 for (x=0; x < (ssize_t) number_grays; x++) 1750 { 1751 /* 1752 Difference variance. 1753 */ 1754 variance.direction[i].red+=density_xy[x].direction[i].red; 1755 variance.direction[i].green+=density_xy[x].direction[i].green; 1756 variance.direction[i].blue+=density_xy[x].direction[i].blue; 1757 if (image->colorspace == CMYKColorspace) 1758 variance.direction[i].black+=density_xy[x].direction[i].black; 1759 if (image->alpha_trait == BlendPixelTrait) 1760 variance.direction[i].alpha+=density_xy[x].direction[i].alpha; 1761 sum_squares.direction[i].red+=density_xy[x].direction[i].red* 1762 density_xy[x].direction[i].red; 1763 sum_squares.direction[i].green+=density_xy[x].direction[i].green* 1764 density_xy[x].direction[i].green; 1765 sum_squares.direction[i].blue+=density_xy[x].direction[i].blue* 1766 density_xy[x].direction[i].blue; 1767 if (image->colorspace == CMYKColorspace) 1768 sum_squares.direction[i].black+=density_xy[x].direction[i].black* 1769 density_xy[x].direction[i].black; 1770 if (image->alpha_trait == BlendPixelTrait) 1771 sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha* 1772 density_xy[x].direction[i].alpha; 1773 /* 1774 Difference entropy. 1775 */ 1776 channel_features[RedPixelChannel].difference_entropy[i]-= 1777 density_xy[x].direction[i].red* 1778 MagickLog10(density_xy[x].direction[i].red); 1779 channel_features[GreenPixelChannel].difference_entropy[i]-= 1780 density_xy[x].direction[i].green* 1781 MagickLog10(density_xy[x].direction[i].green); 1782 channel_features[BluePixelChannel].difference_entropy[i]-= 1783 density_xy[x].direction[i].blue* 1784 MagickLog10(density_xy[x].direction[i].blue); 1785 if (image->colorspace == CMYKColorspace) 1786 channel_features[BlackPixelChannel].difference_entropy[i]-= 1787 density_xy[x].direction[i].black* 1788 MagickLog10(density_xy[x].direction[i].black); 1789 if (image->alpha_trait == BlendPixelTrait) 1790 channel_features[AlphaPixelChannel].difference_entropy[i]-= 1791 density_xy[x].direction[i].alpha* 1792 MagickLog10(density_xy[x].direction[i].alpha); 1793 /* 1794 Information Measures of Correlation. 1795 */ 1796 entropy_x.direction[i].red-=(density_x[x].direction[i].red* 1797 MagickLog10(density_x[x].direction[i].red)); 1798 entropy_x.direction[i].green-=(density_x[x].direction[i].green* 1799 MagickLog10(density_x[x].direction[i].green)); 1800 entropy_x.direction[i].blue-=(density_x[x].direction[i].blue* 1801 MagickLog10(density_x[x].direction[i].blue)); 1802 if (image->colorspace == CMYKColorspace) 1803 entropy_x.direction[i].black-=(density_x[x].direction[i].black* 1804 MagickLog10(density_x[x].direction[i].black)); 1805 if (image->alpha_trait == BlendPixelTrait) 1806 entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha* 1807 MagickLog10(density_x[x].direction[i].alpha)); 1808 entropy_y.direction[i].red-=(density_y[x].direction[i].red* 1809 MagickLog10(density_y[x].direction[i].red)); 1810 entropy_y.direction[i].green-=(density_y[x].direction[i].green* 1811 MagickLog10(density_y[x].direction[i].green)); 1812 entropy_y.direction[i].blue-=(density_y[x].direction[i].blue* 1813 MagickLog10(density_y[x].direction[i].blue)); 1814 if (image->colorspace == CMYKColorspace) 1815 entropy_y.direction[i].black-=(density_y[x].direction[i].black* 1816 MagickLog10(density_y[x].direction[i].black)); 1817 if (image->alpha_trait == BlendPixelTrait) 1818 entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha* 1819 MagickLog10(density_y[x].direction[i].alpha)); 1820 } 1821 /* 1822 Difference variance. 1823 */ 1824 channel_features[RedPixelChannel].difference_variance[i]= 1825 (((double) number_grays*number_grays*sum_squares.direction[i].red)- 1826 (variance.direction[i].red*variance.direction[i].red))/ 1827 ((double) number_grays*number_grays*number_grays*number_grays); 1828 channel_features[GreenPixelChannel].difference_variance[i]= 1829 (((double) number_grays*number_grays*sum_squares.direction[i].green)- 1830 (variance.direction[i].green*variance.direction[i].green))/ 1831 ((double) number_grays*number_grays*number_grays*number_grays); 1832 channel_features[BluePixelChannel].difference_variance[i]= 1833 (((double) number_grays*number_grays*sum_squares.direction[i].blue)- 1834 (variance.direction[i].blue*variance.direction[i].blue))/ 1835 ((double) number_grays*number_grays*number_grays*number_grays); 1836 if (image->colorspace == CMYKColorspace) 1837 channel_features[BlackPixelChannel].difference_variance[i]= 1838 (((double) number_grays*number_grays*sum_squares.direction[i].black)- 1839 (variance.direction[i].black*variance.direction[i].black))/ 1840 ((double) number_grays*number_grays*number_grays*number_grays); 1841 if (image->alpha_trait == BlendPixelTrait) 1842 channel_features[AlphaPixelChannel].difference_variance[i]= 1843 (((double) number_grays*number_grays*sum_squares.direction[i].alpha)- 1844 (variance.direction[i].alpha*variance.direction[i].alpha))/ 1845 ((double) number_grays*number_grays*number_grays*number_grays); 1846 /* 1847 Information Measures of Correlation. 1848 */ 1849 channel_features[RedPixelChannel].measure_of_correlation_1[i]= 1850 (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/ 1851 (entropy_x.direction[i].red > entropy_y.direction[i].red ? 1852 entropy_x.direction[i].red : entropy_y.direction[i].red); 1853 channel_features[GreenPixelChannel].measure_of_correlation_1[i]= 1854 (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/ 1855 (entropy_x.direction[i].green > entropy_y.direction[i].green ? 1856 entropy_x.direction[i].green : entropy_y.direction[i].green); 1857 channel_features[BluePixelChannel].measure_of_correlation_1[i]= 1858 (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/ 1859 (entropy_x.direction[i].blue > entropy_y.direction[i].blue ? 1860 entropy_x.direction[i].blue : entropy_y.direction[i].blue); 1861 if (image->colorspace == CMYKColorspace) 1862 channel_features[BlackPixelChannel].measure_of_correlation_1[i]= 1863 (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/ 1864 (entropy_x.direction[i].black > entropy_y.direction[i].black ? 1865 entropy_x.direction[i].black : entropy_y.direction[i].black); 1866 if (image->alpha_trait == BlendPixelTrait) 1867 channel_features[AlphaPixelChannel].measure_of_correlation_1[i]= 1868 (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/ 1869 (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ? 1870 entropy_x.direction[i].alpha : entropy_y.direction[i].alpha); 1871 channel_features[RedPixelChannel].measure_of_correlation_2[i]= 1872 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red- 1873 entropy_xy.direction[i].red))))); 1874 channel_features[GreenPixelChannel].measure_of_correlation_2[i]= 1875 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green- 1876 entropy_xy.direction[i].green))))); 1877 channel_features[BluePixelChannel].measure_of_correlation_2[i]= 1878 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue- 1879 entropy_xy.direction[i].blue))))); 1880 if (image->colorspace == CMYKColorspace) 1881 channel_features[BlackPixelChannel].measure_of_correlation_2[i]= 1882 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black- 1883 entropy_xy.direction[i].black))))); 1884 if (image->alpha_trait == BlendPixelTrait) 1885 channel_features[AlphaPixelChannel].measure_of_correlation_2[i]= 1886 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha- 1887 entropy_xy.direction[i].alpha))))); 1888 } 1889 /* 1890 Compute more texture features. 1891 */ 1892#if defined(MAGICKCORE_OPENMP_SUPPORT) 1893 #pragma omp parallel for schedule(static,4) shared(status) \ 1894 magick_threads(image,image,number_grays,1) 1895#endif 1896 for (i=0; i < 4; i++) 1897 { 1898 ssize_t 1899 z; 1900 1901 for (z=0; z < (ssize_t) number_grays; z++) 1902 { 1903 register ssize_t 1904 y; 1905 1906 ChannelStatistics 1907 pixel; 1908 1909 (void) ResetMagickMemory(&pixel,0,sizeof(pixel)); 1910 for (y=0; y < (ssize_t) number_grays; y++) 1911 { 1912 register ssize_t 1913 x; 1914 1915 for (x=0; x < (ssize_t) number_grays; x++) 1916 { 1917 /* 1918 Contrast: amount of local variations present in an image. 1919 */ 1920 if (((y-x) == z) || ((x-y) == z)) 1921 { 1922 pixel.direction[i].red+=cooccurrence[x][y].direction[i].red; 1923 pixel.direction[i].green+=cooccurrence[x][y].direction[i].green; 1924 pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue; 1925 if (image->colorspace == CMYKColorspace) 1926 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black; 1927 if (image->alpha_trait == BlendPixelTrait) 1928 pixel.direction[i].alpha+= 1929 cooccurrence[x][y].direction[i].alpha; 1930 } 1931 /* 1932 Maximum Correlation Coefficient. 1933 */ 1934 Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red* 1935 cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/ 1936 density_y[x].direction[i].red; 1937 Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green* 1938 cooccurrence[y][x].direction[i].green/ 1939 density_x[z].direction[i].green/density_y[x].direction[i].red; 1940 Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue* 1941 cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/ 1942 density_y[x].direction[i].blue; 1943 if (image->colorspace == CMYKColorspace) 1944 Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black* 1945 cooccurrence[y][x].direction[i].black/ 1946 density_x[z].direction[i].black/density_y[x].direction[i].black; 1947 if (image->alpha_trait == BlendPixelTrait) 1948 Q[z][y].direction[i].alpha+= 1949 cooccurrence[z][x].direction[i].alpha* 1950 cooccurrence[y][x].direction[i].alpha/ 1951 density_x[z].direction[i].alpha/ 1952 density_y[x].direction[i].alpha; 1953 } 1954 } 1955 channel_features[RedPixelChannel].contrast[i]+=z*z* 1956 pixel.direction[i].red; 1957 channel_features[GreenPixelChannel].contrast[i]+=z*z* 1958 pixel.direction[i].green; 1959 channel_features[BluePixelChannel].contrast[i]+=z*z* 1960 pixel.direction[i].blue; 1961 if (image->colorspace == CMYKColorspace) 1962 channel_features[BlackPixelChannel].contrast[i]+=z*z* 1963 pixel.direction[i].black; 1964 if (image->alpha_trait == BlendPixelTrait) 1965 channel_features[AlphaPixelChannel].contrast[i]+=z*z* 1966 pixel.direction[i].alpha; 1967 } 1968 /* 1969 Maximum Correlation Coefficient. 1970 Future: return second largest eigenvalue of Q. 1971 */ 1972 channel_features[RedPixelChannel].maximum_correlation_coefficient[i]= 1973 sqrt((double) -1.0); 1974 channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]= 1975 sqrt((double) -1.0); 1976 channel_features[BluePixelChannel].maximum_correlation_coefficient[i]= 1977 sqrt((double) -1.0); 1978 if (image->colorspace == CMYKColorspace) 1979 channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]= 1980 sqrt((double) -1.0); 1981 if (image->alpha_trait == BlendPixelTrait) 1982 channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]= 1983 sqrt((double) -1.0); 1984 } 1985 /* 1986 Relinquish resources. 1987 */ 1988 sum=(ChannelStatistics *) RelinquishMagickMemory(sum); 1989 for (i=0; i < (ssize_t) number_grays; i++) 1990 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]); 1991 Q=(ChannelStatistics **) RelinquishMagickMemory(Q); 1992 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y); 1993 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy); 1994 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x); 1995 for (i=0; i < (ssize_t) number_grays; i++) 1996 cooccurrence[i]=(ChannelStatistics *) 1997 RelinquishMagickMemory(cooccurrence[i]); 1998 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence); 1999 return(channel_features); 2000} 2001