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