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