morphology.c revision 306b0f1c3332736bfaddd94c5aaff078a69cbac3
1/* 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3% % 4% % 5% % 6% M M OOO RRRR PPPP H H OOO L OOO GGGG Y Y % 7% MM MM O O R R P P H H O O L O O G Y Y % 8% M M M O O RRRR PPPP HHHHH O O L O O G GGG Y % 9% M M O O R R P H H O O L O O G G Y % 10% M M OOO R R P H H OOO LLLLL OOO GGG Y % 11% % 12% % 13% MagickCore Morphology Methods % 14% % 15% Software Design % 16% Anthony Thyssen % 17% January 2010 % 18% % 19% % 20% Copyright 1999-2013 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% Morpology is the the application of various kernels, of any size and even 37% shape, to a image in various ways (typically binary, but not always). 38% 39% Convolution (weighted sum or average) is just one specific type of 40% morphology. Just one that is very common for image bluring and sharpening 41% effects. Not only 2D Gaussian blurring, but also 2-pass 1D Blurring. 42% 43% This module provides not only a general morphology function, and the ability 44% to apply more advanced or iterative morphologies, but also functions for the 45% generation of many different types of kernel arrays from user supplied 46% arguments. Prehaps even the generation of a kernel from a small image. 47*/ 48 49/* 50 Include declarations. 51*/ 52#include "MagickCore/studio.h" 53#include "MagickCore/artifact.h" 54#include "MagickCore/cache-view.h" 55#include "MagickCore/color-private.h" 56#include "MagickCore/enhance.h" 57#include "MagickCore/exception.h" 58#include "MagickCore/exception-private.h" 59#include "MagickCore/gem.h" 60#include "MagickCore/gem-private.h" 61#include "MagickCore/hashmap.h" 62#include "MagickCore/image.h" 63#include "MagickCore/image-private.h" 64#include "MagickCore/list.h" 65#include "MagickCore/magick.h" 66#include "MagickCore/memory_.h" 67#include "MagickCore/memory-private.h" 68#include "MagickCore/monitor-private.h" 69#include "MagickCore/morphology.h" 70#include "MagickCore/morphology-private.h" 71#include "MagickCore/option.h" 72#include "MagickCore/pixel-accessor.h" 73#include "MagickCore/pixel-private.h" 74#include "MagickCore/prepress.h" 75#include "MagickCore/quantize.h" 76#include "MagickCore/resource_.h" 77#include "MagickCore/registry.h" 78#include "MagickCore/semaphore.h" 79#include "MagickCore/splay-tree.h" 80#include "MagickCore/statistic.h" 81#include "MagickCore/string_.h" 82#include "MagickCore/string-private.h" 83#include "MagickCore/thread-private.h" 84#include "MagickCore/token.h" 85#include "MagickCore/utility.h" 86#include "MagickCore/utility-private.h" 87 88/* 89 Other global definitions used by module. 90*/ 91static inline double MagickMin(const double x,const double y) 92{ 93 return( x < y ? x : y); 94} 95static inline double MagickMax(const double x,const double y) 96{ 97 return( x > y ? x : y); 98} 99#define Minimize(assign,value) assign=MagickMin(assign,value) 100#define Maximize(assign,value) assign=MagickMax(assign,value) 101 102/* Integer Factorial Function - for a Binomial kernel */ 103#if 1 104static inline size_t fact(size_t n) 105{ 106 size_t f,l; 107 for(f=1, l=2; l <= n; f=f*l, l++); 108 return(f); 109} 110#elif 1 /* glibc floating point alternatives */ 111#define fact(n) ((size_t)tgamma((double)n+1)) 112#else 113#define fact(n) ((size_t)lgamma((double)n+1)) 114#endif 115 116 117/* Currently these are only internal to this module */ 118static void 119 CalcKernelMetaData(KernelInfo *), 120 ExpandMirrorKernelInfo(KernelInfo *), 121 ExpandRotateKernelInfo(KernelInfo *, const double), 122 RotateKernelInfo(KernelInfo *, double); 123 124 125/* Quick function to find last kernel in a kernel list */ 126static inline KernelInfo *LastKernelInfo(KernelInfo *kernel) 127{ 128 while (kernel->next != (KernelInfo *) NULL) 129 kernel = kernel->next; 130 return(kernel); 131} 132 133/* 134%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 135% % 136% % 137% % 138% A c q u i r e K e r n e l I n f o % 139% % 140% % 141% % 142%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 143% 144% AcquireKernelInfo() takes the given string (generally supplied by the 145% user) and converts it into a Morphology/Convolution Kernel. This allows 146% users to specify a kernel from a number of pre-defined kernels, or to fully 147% specify their own kernel for a specific Convolution or Morphology 148% Operation. 149% 150% The kernel so generated can be any rectangular array of floating point 151% values (doubles) with the 'control point' or 'pixel being affected' 152% anywhere within that array of values. 153% 154% Previously IM was restricted to a square of odd size using the exact 155% center as origin, this is no longer the case, and any rectangular kernel 156% with any value being declared the origin. This in turn allows the use of 157% highly asymmetrical kernels. 158% 159% The floating point values in the kernel can also include a special value 160% known as 'nan' or 'not a number' to indicate that this value is not part 161% of the kernel array. This allows you to shaped the kernel within its 162% rectangular area. That is 'nan' values provide a 'mask' for the kernel 163% shape. However at least one non-nan value must be provided for correct 164% working of a kernel. 165% 166% The returned kernel should be freed using the DestroyKernelInfo() when you 167% are finished with it. Do not free this memory yourself. 168% 169% Input kernel defintion strings can consist of any of three types. 170% 171% "name:args[[@><]" 172% Select from one of the built in kernels, using the name and 173% geometry arguments supplied. See AcquireKernelBuiltIn() 174% 175% "WxH[+X+Y][@><]:num, num, num ..." 176% a kernel of size W by H, with W*H floating point numbers following. 177% the 'center' can be optionally be defined at +X+Y (such that +0+0 178% is top left corner). If not defined the pixel in the center, for 179% odd sizes, or to the immediate top or left of center for even sizes 180% is automatically selected. 181% 182% "num, num, num, num, ..." 183% list of floating point numbers defining an 'old style' odd sized 184% square kernel. At least 9 values should be provided for a 3x3 185% square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc. 186% Values can be space or comma separated. This is not recommended. 187% 188% You can define a 'list of kernels' which can be used by some morphology 189% operators A list is defined as a semi-colon separated list kernels. 190% 191% " kernel ; kernel ; kernel ; " 192% 193% Any extra ';' characters, at start, end or between kernel defintions are 194% simply ignored. 195% 196% The special flags will expand a single kernel, into a list of rotated 197% kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree 198% cyclic rotations, while a '>' will generate a list of 90-degree rotations. 199% The '<' also exands using 90-degree rotates, but giving a 180-degree 200% reflected kernel before the +/- 90-degree rotations, which can be important 201% for Thinning operations. 202% 203% Note that 'name' kernels will start with an alphabetic character while the 204% new kernel specification has a ':' character in its specification string. 205% If neither is the case, it is assumed an old style of a simple list of 206% numbers generating a odd-sized square kernel has been given. 207% 208% The format of the AcquireKernal method is: 209% 210% KernelInfo *AcquireKernelInfo(const char *kernel_string) 211% 212% A description of each parameter follows: 213% 214% o kernel_string: the Morphology/Convolution kernel wanted. 215% 216*/ 217 218/* This was separated so that it could be used as a separate 219** array input handling function, such as for -color-matrix 220*/ 221static KernelInfo *ParseKernelArray(const char *kernel_string) 222{ 223 KernelInfo 224 *kernel; 225 226 char 227 token[MaxTextExtent]; 228 229 const char 230 *p, 231 *end; 232 233 register ssize_t 234 i; 235 236 double 237 nan = sqrt((double)-1.0); /* Special Value : Not A Number */ 238 239 MagickStatusType 240 flags; 241 242 GeometryInfo 243 args; 244 245 kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel)); 246 if (kernel == (KernelInfo *)NULL) 247 return(kernel); 248 (void) ResetMagickMemory(kernel,0,sizeof(*kernel)); 249 kernel->minimum = kernel->maximum = kernel->angle = 0.0; 250 kernel->negative_range = kernel->positive_range = 0.0; 251 kernel->type = UserDefinedKernel; 252 kernel->next = (KernelInfo *) NULL; 253 kernel->signature = MagickSignature; 254 if (kernel_string == (const char *) NULL) 255 return(kernel); 256 257 /* find end of this specific kernel definition string */ 258 end = strchr(kernel_string, ';'); 259 if ( end == (char *) NULL ) 260 end = strchr(kernel_string, '\0'); 261 262 /* clear flags - for Expanding kernel lists thorugh rotations */ 263 flags = NoValue; 264 265 /* Has a ':' in argument - New user kernel specification 266 FUTURE: this split on ':' could be done by StringToken() 267 */ 268 p = strchr(kernel_string, ':'); 269 if ( p != (char *) NULL && p < end) 270 { 271 /* ParseGeometry() needs the geometry separated! -- Arrgghh */ 272 memcpy(token, kernel_string, (size_t) (p-kernel_string)); 273 token[p-kernel_string] = '\0'; 274 SetGeometryInfo(&args); 275 flags = ParseGeometry(token, &args); 276 277 /* Size handling and checks of geometry settings */ 278 if ( (flags & WidthValue) == 0 ) /* if no width then */ 279 args.rho = args.sigma; /* then width = height */ 280 if ( args.rho < 1.0 ) /* if width too small */ 281 args.rho = 1.0; /* then width = 1 */ 282 if ( args.sigma < 1.0 ) /* if height too small */ 283 args.sigma = args.rho; /* then height = width */ 284 kernel->width = (size_t)args.rho; 285 kernel->height = (size_t)args.sigma; 286 287 /* Offset Handling and Checks */ 288 if ( args.xi < 0.0 || args.psi < 0.0 ) 289 return(DestroyKernelInfo(kernel)); 290 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi 291 : (ssize_t) (kernel->width-1)/2; 292 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi 293 : (ssize_t) (kernel->height-1)/2; 294 if ( kernel->x >= (ssize_t) kernel->width || 295 kernel->y >= (ssize_t) kernel->height ) 296 return(DestroyKernelInfo(kernel)); 297 298 p++; /* advance beyond the ':' */ 299 } 300 else 301 { /* ELSE - Old old specification, forming odd-square kernel */ 302 /* count up number of values given */ 303 p=(const char *) kernel_string; 304 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\'')) 305 p++; /* ignore "'" chars for convolve filter usage - Cristy */ 306 for (i=0; p < end; i++) 307 { 308 GetMagickToken(p,&p,token); 309 if (*token == ',') 310 GetMagickToken(p,&p,token); 311 } 312 /* set the size of the kernel - old sized square */ 313 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0); 314 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 315 p=(const char *) kernel_string; 316 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\'')) 317 p++; /* ignore "'" chars for convolve filter usage - Cristy */ 318 } 319 320 /* Read in the kernel values from rest of input string argument */ 321 kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory( 322 kernel->width,kernel->height*sizeof(*kernel->values))); 323 if (kernel->values == (MagickRealType *) NULL) 324 return(DestroyKernelInfo(kernel)); 325 kernel->minimum = +MagickHuge; 326 kernel->maximum = -MagickHuge; 327 kernel->negative_range = kernel->positive_range = 0.0; 328 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++) 329 { 330 GetMagickToken(p,&p,token); 331 if (*token == ',') 332 GetMagickToken(p,&p,token); 333 if ( LocaleCompare("nan",token) == 0 334 || LocaleCompare("-",token) == 0 ) { 335 kernel->values[i] = nan; /* this value is not part of neighbourhood */ 336 } 337 else { 338 kernel->values[i] = StringToDouble(token,(char **) NULL); 339 ( kernel->values[i] < 0) 340 ? ( kernel->negative_range += kernel->values[i] ) 341 : ( kernel->positive_range += kernel->values[i] ); 342 Minimize(kernel->minimum, kernel->values[i]); 343 Maximize(kernel->maximum, kernel->values[i]); 344 } 345 } 346 347 /* sanity check -- no more values in kernel definition */ 348 GetMagickToken(p,&p,token); 349 if ( *token != '\0' && *token != ';' && *token != '\'' ) 350 return(DestroyKernelInfo(kernel)); 351 352#if 0 353 /* this was the old method of handling a incomplete kernel */ 354 if ( i < (ssize_t) (kernel->width*kernel->height) ) { 355 Minimize(kernel->minimum, kernel->values[i]); 356 Maximize(kernel->maximum, kernel->values[i]); 357 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++) 358 kernel->values[i]=0.0; 359 } 360#else 361 /* Number of values for kernel was not enough - Report Error */ 362 if ( i < (ssize_t) (kernel->width*kernel->height) ) 363 return(DestroyKernelInfo(kernel)); 364#endif 365 366 /* check that we recieved at least one real (non-nan) value! */ 367 if ( kernel->minimum == MagickHuge ) 368 return(DestroyKernelInfo(kernel)); 369 370 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */ 371 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */ 372 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */ 373 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */ 374 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */ 375 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */ 376 377 return(kernel); 378} 379 380static KernelInfo *ParseKernelName(const char *kernel_string) 381{ 382 char 383 token[MaxTextExtent]; 384 385 const char 386 *p, 387 *end; 388 389 GeometryInfo 390 args; 391 392 KernelInfo 393 *kernel; 394 395 MagickStatusType 396 flags; 397 398 ssize_t 399 type; 400 401 /* Parse special 'named' kernel */ 402 GetMagickToken(kernel_string,&p,token); 403 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token); 404 if ( type < 0 || type == UserDefinedKernel ) 405 return((KernelInfo *)NULL); /* not a valid named kernel */ 406 407 while (((isspace((int) ((unsigned char) *p)) != 0) || 408 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';')) 409 p++; 410 411 end = strchr(p, ';'); /* end of this kernel defintion */ 412 if ( end == (char *) NULL ) 413 end = strchr(p, '\0'); 414 415 /* ParseGeometry() needs the geometry separated! -- Arrgghh */ 416 memcpy(token, p, (size_t) (end-p)); 417 token[end-p] = '\0'; 418 SetGeometryInfo(&args); 419 flags = ParseGeometry(token, &args); 420 421#if 0 422 /* For Debugging Geometry Input */ 423 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n", 424 flags, args.rho, args.sigma, args.xi, args.psi ); 425#endif 426 427 /* special handling of missing values in input string */ 428 switch( type ) { 429 /* Shape Kernel Defaults */ 430 case UnityKernel: 431 if ( (flags & WidthValue) == 0 ) 432 args.rho = 1.0; /* Default scale = 1.0, zero is valid */ 433 break; 434 case SquareKernel: 435 case DiamondKernel: 436 case OctagonKernel: 437 case DiskKernel: 438 case PlusKernel: 439 case CrossKernel: 440 if ( (flags & HeightValue) == 0 ) 441 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */ 442 break; 443 case RingKernel: 444 if ( (flags & XValue) == 0 ) 445 args.xi = 1.0; /* Default scale = 1.0, zero is valid */ 446 break; 447 case RectangleKernel: /* Rectangle - set size defaults */ 448 if ( (flags & WidthValue) == 0 ) /* if no width then */ 449 args.rho = args.sigma; /* then width = height */ 450 if ( args.rho < 1.0 ) /* if width too small */ 451 args.rho = 3; /* then width = 3 */ 452 if ( args.sigma < 1.0 ) /* if height too small */ 453 args.sigma = args.rho; /* then height = width */ 454 if ( (flags & XValue) == 0 ) /* center offset if not defined */ 455 args.xi = (double)(((ssize_t)args.rho-1)/2); 456 if ( (flags & YValue) == 0 ) 457 args.psi = (double)(((ssize_t)args.sigma-1)/2); 458 break; 459 /* Distance Kernel Defaults */ 460 case ChebyshevKernel: 461 case ManhattanKernel: 462 case OctagonalKernel: 463 case EuclideanKernel: 464 if ( (flags & HeightValue) == 0 ) /* no distance scale */ 465 args.sigma = 100.0; /* default distance scaling */ 466 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */ 467 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */ 468 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */ 469 args.sigma *= QuantumRange/100.0; /* percentage of color range */ 470 break; 471 default: 472 break; 473 } 474 475 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args); 476 if ( kernel == (KernelInfo *) NULL ) 477 return(kernel); 478 479 /* global expand to rotated kernel list - only for single kernels */ 480 if ( kernel->next == (KernelInfo *) NULL ) { 481 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */ 482 ExpandRotateKernelInfo(kernel, 45.0); 483 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */ 484 ExpandRotateKernelInfo(kernel, 90.0); 485 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */ 486 ExpandMirrorKernelInfo(kernel); 487 } 488 489 return(kernel); 490} 491 492MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string) 493{ 494 495 KernelInfo 496 *kernel, 497 *new_kernel; 498 499 char 500 token[MaxTextExtent]; 501 502 const char 503 *p; 504 505 size_t 506 kernel_number; 507 508 if (kernel_string == (const char *) NULL) 509 return(ParseKernelArray(kernel_string)); 510 p = kernel_string; 511 kernel = NULL; 512 kernel_number = 0; 513 514 while ( GetMagickToken(p,NULL,token), *token != '\0' ) { 515 516 /* ignore extra or multiple ';' kernel separators */ 517 if ( *token != ';' ) { 518 519 /* tokens starting with alpha is a Named kernel */ 520 if (isalpha((int) *token) != 0) 521 new_kernel = ParseKernelName(p); 522 else /* otherwise a user defined kernel array */ 523 new_kernel = ParseKernelArray(p); 524 525 /* Error handling -- this is not proper error handling! */ 526 if ( new_kernel == (KernelInfo *) NULL ) { 527 (void) FormatLocaleFile(stderr,"Failed to parse kernel number #%.20g\n", 528 (double) kernel_number); 529 if ( kernel != (KernelInfo *) NULL ) 530 kernel=DestroyKernelInfo(kernel); 531 return((KernelInfo *) NULL); 532 } 533 534 /* initialise or append the kernel list */ 535 if ( kernel == (KernelInfo *) NULL ) 536 kernel = new_kernel; 537 else 538 LastKernelInfo(kernel)->next = new_kernel; 539 } 540 541 /* look for the next kernel in list */ 542 p = strchr(p, ';'); 543 if ( p == (char *) NULL ) 544 break; 545 p++; 546 547 } 548 return(kernel); 549} 550 551 552/* 553%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 554% % 555% % 556% % 557% A c q u i r e K e r n e l B u i l t I n % 558% % 559% % 560% % 561%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 562% 563% AcquireKernelBuiltIn() returned one of the 'named' built-in types of 564% kernels used for special purposes such as gaussian blurring, skeleton 565% pruning, and edge distance determination. 566% 567% They take a KernelType, and a set of geometry style arguments, which were 568% typically decoded from a user supplied string, or from a more complex 569% Morphology Method that was requested. 570% 571% The format of the AcquireKernalBuiltIn method is: 572% 573% KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type, 574% const GeometryInfo args) 575% 576% A description of each parameter follows: 577% 578% o type: the pre-defined type of kernel wanted 579% 580% o args: arguments defining or modifying the kernel 581% 582% Convolution Kernels 583% 584% Unity 585% The a No-Op or Scaling single element kernel. 586% 587% Gaussian:{radius},{sigma} 588% Generate a two-dimensional gaussian kernel, as used by -gaussian. 589% The sigma for the curve is required. The resulting kernel is 590% normalized, 591% 592% If 'sigma' is zero, you get a single pixel on a field of zeros. 593% 594% NOTE: that the 'radius' is optional, but if provided can limit (clip) 595% the final size of the resulting kernel to a square 2*radius+1 in size. 596% The radius should be at least 2 times that of the sigma value, or 597% sever clipping and aliasing may result. If not given or set to 0 the 598% radius will be determined so as to produce the best minimal error 599% result, which is usally much larger than is normally needed. 600% 601% LoG:{radius},{sigma} 602% "Laplacian of a Gaussian" or "Mexician Hat" Kernel. 603% The supposed ideal edge detection, zero-summing kernel. 604% 605% An alturnative to this kernel is to use a "DoG" with a sigma ratio of 606% approx 1.6 (according to wikipedia). 607% 608% DoG:{radius},{sigma1},{sigma2} 609% "Difference of Gaussians" Kernel. 610% As "Gaussian" but with a gaussian produced by 'sigma2' subtracted 611% from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1. 612% The result is a zero-summing kernel. 613% 614% Blur:{radius},{sigma}[,{angle}] 615% Generates a 1 dimensional or linear gaussian blur, at the angle given 616% (current restricted to orthogonal angles). If a 'radius' is given the 617% kernel is clipped to a width of 2*radius+1. Kernel can be rotated 618% by a 90 degree angle. 619% 620% If 'sigma' is zero, you get a single pixel on a field of zeros. 621% 622% Note that two convolutions with two "Blur" kernels perpendicular to 623% each other, is equivalent to a far larger "Gaussian" kernel with the 624% same sigma value, However it is much faster to apply. This is how the 625% "-blur" operator actually works. 626% 627% Comet:{width},{sigma},{angle} 628% Blur in one direction only, much like how a bright object leaves 629% a comet like trail. The Kernel is actually half a gaussian curve, 630% Adding two such blurs in opposite directions produces a Blur Kernel. 631% Angle can be rotated in multiples of 90 degrees. 632% 633% Note that the first argument is the width of the kernel and not the 634% radius of the kernel. 635% 636% Binomial:[{radius}] 637% Generate a discrete kernel using a 2 dimentional Pascel's Triangle 638% of values. Used for special forma of image filters. 639% 640% # Still to be implemented... 641% # 642% # Filter2D 643% # Filter1D 644% # Set kernel values using a resize filter, and given scale (sigma) 645% # Cylindrical or Linear. Is this possible with an image? 646% # 647% 648% Named Constant Convolution Kernels 649% 650% All these are unscaled, zero-summing kernels by default. As such for 651% non-HDRI version of ImageMagick some form of normalization, user scaling, 652% and biasing the results is recommended, to prevent the resulting image 653% being 'clipped'. 654% 655% The 3x3 kernels (most of these) can be circularly rotated in multiples of 656% 45 degrees to generate the 8 angled varients of each of the kernels. 657% 658% Laplacian:{type} 659% Discrete Lapacian Kernels, (without normalization) 660% Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood) 661% Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood) 662% Type 2 : 3x3 with center:4 edge:1 corner:-2 663% Type 3 : 3x3 with center:4 edge:-2 corner:1 664% Type 5 : 5x5 laplacian 665% Type 7 : 7x7 laplacian 666% Type 15 : 5x5 LoG (sigma approx 1.4) 667% Type 19 : 9x9 LoG (sigma approx 1.4) 668% 669% Sobel:{angle} 670% Sobel 'Edge' convolution kernel (3x3) 671% | -1, 0, 1 | 672% | -2, 0,-2 | 673% | -1, 0, 1 | 674% 675% Roberts:{angle} 676% Roberts convolution kernel (3x3) 677% | 0, 0, 0 | 678% | -1, 1, 0 | 679% | 0, 0, 0 | 680% 681% Prewitt:{angle} 682% Prewitt Edge convolution kernel (3x3) 683% | -1, 0, 1 | 684% | -1, 0, 1 | 685% | -1, 0, 1 | 686% 687% Compass:{angle} 688% Prewitt's "Compass" convolution kernel (3x3) 689% | -1, 1, 1 | 690% | -1,-2, 1 | 691% | -1, 1, 1 | 692% 693% Kirsch:{angle} 694% Kirsch's "Compass" convolution kernel (3x3) 695% | -3,-3, 5 | 696% | -3, 0, 5 | 697% | -3,-3, 5 | 698% 699% FreiChen:{angle} 700% Frei-Chen Edge Detector is based on a kernel that is similar to 701% the Sobel Kernel, but is designed to be isotropic. That is it takes 702% into account the distance of the diagonal in the kernel. 703% 704% | 1, 0, -1 | 705% | sqrt(2), 0, -sqrt(2) | 706% | 1, 0, -1 | 707% 708% FreiChen:{type},{angle} 709% 710% Frei-Chen Pre-weighted kernels... 711% 712% Type 0: default un-nomalized version shown above. 713% 714% Type 1: Orthogonal Kernel (same as type 11 below) 715% | 1, 0, -1 | 716% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2) 717% | 1, 0, -1 | 718% 719% Type 2: Diagonal form of Kernel... 720% | 1, sqrt(2), 0 | 721% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2) 722% | 0, -sqrt(2) -1 | 723% 724% However this kernel is als at the heart of the FreiChen Edge Detection 725% Process which uses a set of 9 specially weighted kernel. These 9 726% kernels not be normalized, but directly applied to the image. The 727% results is then added together, to produce the intensity of an edge in 728% a specific direction. The square root of the pixel value can then be 729% taken as the cosine of the edge, and at least 2 such runs at 90 degrees 730% from each other, both the direction and the strength of the edge can be 731% determined. 732% 733% Type 10: All 9 of the following pre-weighted kernels... 734% 735% Type 11: | 1, 0, -1 | 736% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2) 737% | 1, 0, -1 | 738% 739% Type 12: | 1, sqrt(2), 1 | 740% | 0, 0, 0 | / 2*sqrt(2) 741% | 1, sqrt(2), 1 | 742% 743% Type 13: | sqrt(2), -1, 0 | 744% | -1, 0, 1 | / 2*sqrt(2) 745% | 0, 1, -sqrt(2) | 746% 747% Type 14: | 0, 1, -sqrt(2) | 748% | -1, 0, 1 | / 2*sqrt(2) 749% | sqrt(2), -1, 0 | 750% 751% Type 15: | 0, -1, 0 | 752% | 1, 0, 1 | / 2 753% | 0, -1, 0 | 754% 755% Type 16: | 1, 0, -1 | 756% | 0, 0, 0 | / 2 757% | -1, 0, 1 | 758% 759% Type 17: | 1, -2, 1 | 760% | -2, 4, -2 | / 6 761% | -1, -2, 1 | 762% 763% Type 18: | -2, 1, -2 | 764% | 1, 4, 1 | / 6 765% | -2, 1, -2 | 766% 767% Type 19: | 1, 1, 1 | 768% | 1, 1, 1 | / 3 769% | 1, 1, 1 | 770% 771% The first 4 are for edge detection, the next 4 are for line detection 772% and the last is to add a average component to the results. 773% 774% Using a special type of '-1' will return all 9 pre-weighted kernels 775% as a multi-kernel list, so that you can use them directly (without 776% normalization) with the special "-set option:morphology:compose Plus" 777% setting to apply the full FreiChen Edge Detection Technique. 778% 779% If 'type' is large it will be taken to be an actual rotation angle for 780% the default FreiChen (type 0) kernel. As such FreiChen:45 will look 781% like a Sobel:45 but with 'sqrt(2)' instead of '2' values. 782% 783% WARNING: The above was layed out as per 784% http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf 785% But rotated 90 degrees so direction is from left rather than the top. 786% I have yet to find any secondary confirmation of the above. The only 787% other source found was actual source code at 788% http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf 789% Neigher paper defineds the kernels in a way that looks locical or 790% correct when taken as a whole. 791% 792% Boolean Kernels 793% 794% Diamond:[{radius}[,{scale}]] 795% Generate a diamond shaped kernel with given radius to the points. 796% Kernel size will again be radius*2+1 square and defaults to radius 1, 797% generating a 3x3 kernel that is slightly larger than a square. 798% 799% Square:[{radius}[,{scale}]] 800% Generate a square shaped kernel of size radius*2+1, and defaulting 801% to a 3x3 (radius 1). 802% 803% Octagon:[{radius}[,{scale}]] 804% Generate octagonal shaped kernel of given radius and constant scale. 805% Default radius is 3 producing a 7x7 kernel. A radius of 1 will result 806% in "Diamond" kernel. 807% 808% Disk:[{radius}[,{scale}]] 809% Generate a binary disk, thresholded at the radius given, the radius 810% may be a float-point value. Final Kernel size is floor(radius)*2+1 811% square. A radius of 5.3 is the default. 812% 813% NOTE: That a low radii Disk kernels produce the same results as 814% many of the previously defined kernels, but differ greatly at larger 815% radii. Here is a table of equivalences... 816% "Disk:1" => "Diamond", "Octagon:1", or "Cross:1" 817% "Disk:1.5" => "Square" 818% "Disk:2" => "Diamond:2" 819% "Disk:2.5" => "Octagon" 820% "Disk:2.9" => "Square:2" 821% "Disk:3.5" => "Octagon:3" 822% "Disk:4.5" => "Octagon:4" 823% "Disk:5.4" => "Octagon:5" 824% "Disk:6.4" => "Octagon:6" 825% All other Disk shapes are unique to this kernel, but because a "Disk" 826% is more circular when using a larger radius, using a larger radius is 827% preferred over iterating the morphological operation. 828% 829% Rectangle:{geometry} 830% Simply generate a rectangle of 1's with the size given. You can also 831% specify the location of the 'control point', otherwise the closest 832% pixel to the center of the rectangle is selected. 833% 834% Properly centered and odd sized rectangles work the best. 835% 836% Symbol Dilation Kernels 837% 838% These kernel is not a good general morphological kernel, but is used 839% more for highlighting and marking any single pixels in an image using, 840% a "Dilate" method as appropriate. 841% 842% For the same reasons iterating these kernels does not produce the 843% same result as using a larger radius for the symbol. 844% 845% Plus:[{radius}[,{scale}]] 846% Cross:[{radius}[,{scale}]] 847% Generate a kernel in the shape of a 'plus' or a 'cross' with 848% a each arm the length of the given radius (default 2). 849% 850% NOTE: "plus:1" is equivalent to a "Diamond" kernel. 851% 852% Ring:{radius1},{radius2}[,{scale}] 853% A ring of the values given that falls between the two radii. 854% Defaults to a ring of approximataly 3 radius in a 7x7 kernel. 855% This is the 'edge' pixels of the default "Disk" kernel, 856% More specifically, "Ring" -> "Ring:2.5,3.5,1.0" 857% 858% Hit and Miss Kernels 859% 860% Peak:radius1,radius2 861% Find any peak larger than the pixels the fall between the two radii. 862% The default ring of pixels is as per "Ring". 863% Edges 864% Find flat orthogonal edges of a binary shape 865% Corners 866% Find 90 degree corners of a binary shape 867% Diagonals:type 868% A special kernel to thin the 'outside' of diagonals 869% LineEnds:type 870% Find end points of lines (for pruning a skeletion) 871% Two types of lines ends (default to both) can be searched for 872% Type 0: All line ends 873% Type 1: single kernel for 4-conneected line ends 874% Type 2: single kernel for simple line ends 875% LineJunctions 876% Find three line junctions (within a skeletion) 877% Type 0: all line junctions 878% Type 1: Y Junction kernel 879% Type 2: Diagonal T Junction kernel 880% Type 3: Orthogonal T Junction kernel 881% Type 4: Diagonal X Junction kernel 882% Type 5: Orthogonal + Junction kernel 883% Ridges:type 884% Find single pixel ridges or thin lines 885% Type 1: Fine single pixel thick lines and ridges 886% Type 2: Find two pixel thick lines and ridges 887% ConvexHull 888% Octagonal Thickening Kernel, to generate convex hulls of 45 degrees 889% Skeleton:type 890% Traditional skeleton generating kernels. 891% Type 1: Tradional Skeleton kernel (4 connected skeleton) 892% Type 2: HIPR2 Skeleton kernel (8 connected skeleton) 893% Type 3: Thinning skeleton based on a ressearch paper by 894% Dan S. Bloomberg (Default Type) 895% ThinSE:type 896% A huge variety of Thinning Kernels designed to preserve conectivity. 897% many other kernel sets use these kernels as source definitions. 898% Type numbers are 41-49, 81-89, 481, and 482 which are based on 899% the super and sub notations used in the source research paper. 900% 901% Distance Measuring Kernels 902% 903% Different types of distance measuring methods, which are used with the 904% a 'Distance' morphology method for generating a gradient based on 905% distance from an edge of a binary shape, though there is a technique 906% for handling a anti-aliased shape. 907% 908% See the 'Distance' Morphological Method, for information of how it is 909% applied. 910% 911% Chebyshev:[{radius}][x{scale}[%!]] 912% Chebyshev Distance (also known as Tchebychev or Chessboard distance) 913% is a value of one to any neighbour, orthogonal or diagonal. One why 914% of thinking of it is the number of squares a 'King' or 'Queen' in 915% chess needs to traverse reach any other position on a chess board. 916% It results in a 'square' like distance function, but one where 917% diagonals are given a value that is closer than expected. 918% 919% Manhattan:[{radius}][x{scale}[%!]] 920% Manhattan Distance (also known as Rectilinear, City Block, or the Taxi 921% Cab distance metric), it is the distance needed when you can only 922% travel in horizontal or vertical directions only. It is the 923% distance a 'Rook' in chess would have to travel, and results in a 924% diamond like distances, where diagonals are further than expected. 925% 926% Octagonal:[{radius}][x{scale}[%!]] 927% An interleving of Manhatten and Chebyshev metrics producing an 928% increasing octagonally shaped distance. Distances matches those of 929% the "Octagon" shaped kernel of the same radius. The minimum radius 930% and default is 2, producing a 5x5 kernel. 931% 932% Euclidean:[{radius}][x{scale}[%!]] 933% Euclidean distance is the 'direct' or 'as the crow flys' distance. 934% However by default the kernel size only has a radius of 1, which 935% limits the distance to 'Knight' like moves, with only orthogonal and 936% diagonal measurements being correct. As such for the default kernel 937% you will get octagonal like distance function. 938% 939% However using a larger radius such as "Euclidean:4" you will get a 940% much smoother distance gradient from the edge of the shape. Especially 941% if the image is pre-processed to include any anti-aliasing pixels. 942% Of course a larger kernel is slower to use, and not always needed. 943% 944% The first three Distance Measuring Kernels will only generate distances 945% of exact multiples of {scale} in binary images. As such you can use a 946% scale of 1 without loosing any information. However you also need some 947% scaling when handling non-binary anti-aliased shapes. 948% 949% The "Euclidean" Distance Kernel however does generate a non-integer 950% fractional results, and as such scaling is vital even for binary shapes. 951% 952*/ 953 954MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type, 955 const GeometryInfo *args) 956{ 957 KernelInfo 958 *kernel; 959 960 register ssize_t 961 i; 962 963 register ssize_t 964 u, 965 v; 966 967 double 968 nan = sqrt((double)-1.0); /* Special Value : Not A Number */ 969 970 /* Generate a new empty kernel if needed */ 971 kernel=(KernelInfo *) NULL; 972 switch(type) { 973 case UndefinedKernel: /* These should not call this function */ 974 case UserDefinedKernel: 975 assert("Should not call this function" != (char *)NULL); 976 break; 977 case LaplacianKernel: /* Named Descrete Convolution Kernels */ 978 case SobelKernel: /* these are defined using other kernels */ 979 case RobertsKernel: 980 case PrewittKernel: 981 case CompassKernel: 982 case KirschKernel: 983 case FreiChenKernel: 984 case EdgesKernel: /* Hit and Miss kernels */ 985 case CornersKernel: 986 case DiagonalsKernel: 987 case LineEndsKernel: 988 case LineJunctionsKernel: 989 case RidgesKernel: 990 case ConvexHullKernel: 991 case SkeletonKernel: 992 case ThinSEKernel: 993 break; /* A pre-generated kernel is not needed */ 994#if 0 995 /* set to 1 to do a compile-time check that we haven't missed anything */ 996 case UnityKernel: 997 case GaussianKernel: 998 case DoGKernel: 999 case LoGKernel: 1000 case BlurKernel: 1001 case CometKernel: 1002 case BinomialKernel: 1003 case DiamondKernel: 1004 case SquareKernel: 1005 case RectangleKernel: 1006 case OctagonKernel: 1007 case DiskKernel: 1008 case PlusKernel: 1009 case CrossKernel: 1010 case RingKernel: 1011 case PeaksKernel: 1012 case ChebyshevKernel: 1013 case ManhattanKernel: 1014 case OctangonalKernel: 1015 case EuclideanKernel: 1016#else 1017 default: 1018#endif 1019 /* Generate the base Kernel Structure */ 1020 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel)); 1021 if (kernel == (KernelInfo *) NULL) 1022 return(kernel); 1023 (void) ResetMagickMemory(kernel,0,sizeof(*kernel)); 1024 kernel->minimum = kernel->maximum = kernel->angle = 0.0; 1025 kernel->negative_range = kernel->positive_range = 0.0; 1026 kernel->type = type; 1027 kernel->next = (KernelInfo *) NULL; 1028 kernel->signature = MagickSignature; 1029 break; 1030 } 1031 1032 switch(type) { 1033 /* 1034 Convolution Kernels 1035 */ 1036 case UnityKernel: 1037 { 1038 kernel->height = kernel->width = (size_t) 1; 1039 kernel->x = kernel->y = (ssize_t) 0; 1040 kernel->values=(MagickRealType *) MagickAssumeAligned( 1041 AcquireAlignedMemory(1,sizeof(*kernel->values))); 1042 if (kernel->values == (MagickRealType *) NULL) 1043 return(DestroyKernelInfo(kernel)); 1044 kernel->maximum = kernel->values[0] = args->rho; 1045 break; 1046 } 1047 break; 1048 case GaussianKernel: 1049 case DoGKernel: 1050 case LoGKernel: 1051 { double 1052 sigma = fabs(args->sigma), 1053 sigma2 = fabs(args->xi), 1054 A, B, R; 1055 1056 if ( args->rho >= 1.0 ) 1057 kernel->width = (size_t)args->rho*2+1; 1058 else if ( (type != DoGKernel) || (sigma >= sigma2) ) 1059 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma); 1060 else 1061 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2); 1062 kernel->height = kernel->width; 1063 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1064 kernel->values=(MagickRealType *) MagickAssumeAligned( 1065 AcquireAlignedMemory(kernel->width,kernel->height* 1066 sizeof(*kernel->values))); 1067 if (kernel->values == (MagickRealType *) NULL) 1068 return(DestroyKernelInfo(kernel)); 1069 1070 /* WARNING: The following generates a 'sampled gaussian' kernel. 1071 * What we really want is a 'discrete gaussian' kernel. 1072 * 1073 * How to do this is I don't know, but appears to be basied on the 1074 * Error Function 'erf()' (intergral of a gaussian) 1075 */ 1076 1077 if ( type == GaussianKernel || type == DoGKernel ) 1078 { /* Calculate a Gaussian, OR positive half of a DoG */ 1079 if ( sigma > MagickEpsilon ) 1080 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */ 1081 B = (double) (1.0/(Magick2PI*sigma*sigma)); 1082 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 1083 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1084 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B; 1085 } 1086 else /* limiting case - a unity (normalized Dirac) kernel */ 1087 { (void) ResetMagickMemory(kernel->values,0, (size_t) 1088 kernel->width*kernel->height*sizeof(*kernel->values)); 1089 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0; 1090 } 1091 } 1092 1093 if ( type == DoGKernel ) 1094 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */ 1095 if ( sigma2 > MagickEpsilon ) 1096 { sigma = sigma2; /* simplify loop expressions */ 1097 A = 1.0/(2.0*sigma*sigma); 1098 B = (double) (1.0/(Magick2PI*sigma*sigma)); 1099 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 1100 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1101 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B; 1102 } 1103 else /* limiting case - a unity (normalized Dirac) kernel */ 1104 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0; 1105 } 1106 1107 if ( type == LoGKernel ) 1108 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */ 1109 if ( sigma > MagickEpsilon ) 1110 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */ 1111 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma)); 1112 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 1113 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1114 { R = ((double)(u*u+v*v))*A; 1115 kernel->values[i] = (1-R)*exp(-R)*B; 1116 } 1117 } 1118 else /* special case - generate a unity kernel */ 1119 { (void) ResetMagickMemory(kernel->values,0, (size_t) 1120 kernel->width*kernel->height*sizeof(*kernel->values)); 1121 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0; 1122 } 1123 } 1124 1125 /* Note the above kernels may have been 'clipped' by a user defined 1126 ** radius, producing a smaller (darker) kernel. Also for very small 1127 ** sigma's (> 0.1) the central value becomes larger than one, and thus 1128 ** producing a very bright kernel. 1129 ** 1130 ** Normalization will still be needed. 1131 */ 1132 1133 /* Normalize the 2D Gaussian Kernel 1134 ** 1135 ** NB: a CorrelateNormalize performs a normal Normalize if 1136 ** there are no negative values. 1137 */ 1138 CalcKernelMetaData(kernel); /* the other kernel meta-data */ 1139 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue); 1140 1141 break; 1142 } 1143 case BlurKernel: 1144 { double 1145 sigma = fabs(args->sigma), 1146 alpha, beta; 1147 1148 if ( args->rho >= 1.0 ) 1149 kernel->width = (size_t)args->rho*2+1; 1150 else 1151 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma); 1152 kernel->height = 1; 1153 kernel->x = (ssize_t) (kernel->width-1)/2; 1154 kernel->y = 0; 1155 kernel->negative_range = kernel->positive_range = 0.0; 1156 kernel->values=(MagickRealType *) MagickAssumeAligned( 1157 AcquireAlignedMemory(kernel->width,kernel->height* 1158 sizeof(*kernel->values))); 1159 if (kernel->values == (MagickRealType *) NULL) 1160 return(DestroyKernelInfo(kernel)); 1161 1162#if 1 1163#define KernelRank 3 1164 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix). 1165 ** It generates a gaussian 3 times the width, and compresses it into 1166 ** the expected range. This produces a closer normalization of the 1167 ** resulting kernel, especially for very low sigma values. 1168 ** As such while wierd it is prefered. 1169 ** 1170 ** I am told this method originally came from Photoshop. 1171 ** 1172 ** A properly normalized curve is generated (apart from edge clipping) 1173 ** even though we later normalize the result (for edge clipping) 1174 ** to allow the correct generation of a "Difference of Blurs". 1175 */ 1176 1177 /* initialize */ 1178 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */ 1179 (void) ResetMagickMemory(kernel->values,0, (size_t) 1180 kernel->width*kernel->height*sizeof(*kernel->values)); 1181 /* Calculate a Positive 1D Gaussian */ 1182 if ( sigma > MagickEpsilon ) 1183 { sigma *= KernelRank; /* simplify loop expressions */ 1184 alpha = 1.0/(2.0*sigma*sigma); 1185 beta= (double) (1.0/(MagickSQ2PI*sigma )); 1186 for ( u=-v; u <= v; u++) { 1187 kernel->values[(u+v)/KernelRank] += 1188 exp(-((double)(u*u))*alpha)*beta; 1189 } 1190 } 1191 else /* special case - generate a unity kernel */ 1192 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0; 1193#else 1194 /* Direct calculation without curve averaging 1195 This is equivelent to a KernelRank of 1 */ 1196 1197 /* Calculate a Positive Gaussian */ 1198 if ( sigma > MagickEpsilon ) 1199 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */ 1200 beta = 1.0/(MagickSQ2PI*sigma); 1201 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1202 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta; 1203 } 1204 else /* special case - generate a unity kernel */ 1205 { (void) ResetMagickMemory(kernel->values,0, (size_t) 1206 kernel->width*kernel->height*sizeof(*kernel->values)); 1207 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0; 1208 } 1209#endif 1210 /* Note the above kernel may have been 'clipped' by a user defined 1211 ** radius, producing a smaller (darker) kernel. Also for very small 1212 ** sigma's (> 0.1) the central value becomes larger than one, as a 1213 ** result of not generating a actual 'discrete' kernel, and thus 1214 ** producing a very bright 'impulse'. 1215 ** 1216 ** Becuase of these two factors Normalization is required! 1217 */ 1218 1219 /* Normalize the 1D Gaussian Kernel 1220 ** 1221 ** NB: a CorrelateNormalize performs a normal Normalize if 1222 ** there are no negative values. 1223 */ 1224 CalcKernelMetaData(kernel); /* the other kernel meta-data */ 1225 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue); 1226 1227 /* rotate the 1D kernel by given angle */ 1228 RotateKernelInfo(kernel, args->xi ); 1229 break; 1230 } 1231 case CometKernel: 1232 { double 1233 sigma = fabs(args->sigma), 1234 A; 1235 1236 if ( args->rho < 1.0 ) 1237 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1; 1238 else 1239 kernel->width = (size_t)args->rho; 1240 kernel->x = kernel->y = 0; 1241 kernel->height = 1; 1242 kernel->negative_range = kernel->positive_range = 0.0; 1243 kernel->values=(MagickRealType *) MagickAssumeAligned( 1244 AcquireAlignedMemory(kernel->width,kernel->height* 1245 sizeof(*kernel->values))); 1246 if (kernel->values == (MagickRealType *) NULL) 1247 return(DestroyKernelInfo(kernel)); 1248 1249 /* A comet blur is half a 1D gaussian curve, so that the object is 1250 ** blurred in one direction only. This may not be quite the right 1251 ** curve to use so may change in the future. The function must be 1252 ** normalised after generation, which also resolves any clipping. 1253 ** 1254 ** As we are normalizing and not subtracting gaussians, 1255 ** there is no need for a divisor in the gaussian formula 1256 ** 1257 ** It is less comples 1258 */ 1259 if ( sigma > MagickEpsilon ) 1260 { 1261#if 1 1262#define KernelRank 3 1263 v = (ssize_t) kernel->width*KernelRank; /* start/end points */ 1264 (void) ResetMagickMemory(kernel->values,0, (size_t) 1265 kernel->width*sizeof(*kernel->values)); 1266 sigma *= KernelRank; /* simplify the loop expression */ 1267 A = 1.0/(2.0*sigma*sigma); 1268 /* B = 1.0/(MagickSQ2PI*sigma); */ 1269 for ( u=0; u < v; u++) { 1270 kernel->values[u/KernelRank] += 1271 exp(-((double)(u*u))*A); 1272 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */ 1273 } 1274 for (i=0; i < (ssize_t) kernel->width; i++) 1275 kernel->positive_range += kernel->values[i]; 1276#else 1277 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */ 1278 /* B = 1.0/(MagickSQ2PI*sigma); */ 1279 for ( i=0; i < (ssize_t) kernel->width; i++) 1280 kernel->positive_range += 1281 kernel->values[i] = exp(-((double)(i*i))*A); 1282 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */ 1283#endif 1284 } 1285 else /* special case - generate a unity kernel */ 1286 { (void) ResetMagickMemory(kernel->values,0, (size_t) 1287 kernel->width*kernel->height*sizeof(*kernel->values)); 1288 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0; 1289 kernel->positive_range = 1.0; 1290 } 1291 1292 kernel->minimum = 0.0; 1293 kernel->maximum = kernel->values[0]; 1294 kernel->negative_range = 0.0; 1295 1296 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */ 1297 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */ 1298 break; 1299 } 1300 case BinomialKernel: 1301 { 1302 size_t 1303 order_f; 1304 1305 if (args->rho < 1.0) 1306 kernel->width = kernel->height = 3; /* default radius = 1 */ 1307 else 1308 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 1309 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1310 1311 order_f = fact(kernel->width-1); 1312 1313 kernel->values=(MagickRealType *) MagickAssumeAligned( 1314 AcquireAlignedMemory(kernel->width,kernel->height* 1315 sizeof(*kernel->values))); 1316 if (kernel->values == (MagickRealType *) NULL) 1317 return(DestroyKernelInfo(kernel)); 1318 1319 /* set all kernel values within diamond area to scale given */ 1320 for ( i=0, v=0; v < (ssize_t)kernel->height; v++) 1321 { size_t 1322 alpha = order_f / ( fact((size_t) v) * fact(kernel->height-v-1) ); 1323 for ( u=0; u < (ssize_t)kernel->width; u++, i++) 1324 kernel->positive_range += kernel->values[i] = (double) 1325 (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-u-1) )); 1326 } 1327 kernel->minimum = 1.0; 1328 kernel->maximum = kernel->values[kernel->x+kernel->y*kernel->width]; 1329 kernel->negative_range = 0.0; 1330 break; 1331 } 1332 1333 /* 1334 Convolution Kernels - Well Known Named Constant Kernels 1335 */ 1336 case LaplacianKernel: 1337 { switch ( (int) args->rho ) { 1338 case 0: 1339 default: /* laplacian square filter -- default */ 1340 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1"); 1341 break; 1342 case 1: /* laplacian diamond filter */ 1343 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0"); 1344 break; 1345 case 2: 1346 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2"); 1347 break; 1348 case 3: 1349 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1"); 1350 break; 1351 case 5: /* a 5x5 laplacian */ 1352 kernel=ParseKernelArray( 1353 "5: -4,-1,0,-1,-4 -1,2,3,2,-1 0,3,4,3,0 -1,2,3,2,-1 -4,-1,0,-1,-4"); 1354 break; 1355 case 7: /* a 7x7 laplacian */ 1356 kernel=ParseKernelArray( 1357 "7:-10,-5,-2,-1,-2,-5,-10 -5,0,3,4,3,0,-5 -2,3,6,7,6,3,-2 -1,4,7,8,7,4,-1 -2,3,6,7,6,3,-2 -5,0,3,4,3,0,-5 -10,-5,-2,-1,-2,-5,-10" ); 1358 break; 1359 case 15: /* a 5x5 LoG (sigma approx 1.4) */ 1360 kernel=ParseKernelArray( 1361 "5: 0,0,-1,0,0 0,-1,-2,-1,0 -1,-2,16,-2,-1 0,-1,-2,-1,0 0,0,-1,0,0"); 1362 break; 1363 case 19: /* a 9x9 LoG (sigma approx 1.4) */ 1364 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */ 1365 kernel=ParseKernelArray( 1366 "9: 0,-1,-1,-2,-2,-2,-1,-1,0 -1,-2,-4,-5,-5,-5,-4,-2,-1 -1,-4,-5,-3,-0,-3,-5,-4,-1 -2,-5,-3,12,24,12,-3,-5,-2 -2,-5,-0,24,40,24,-0,-5,-2 -2,-5,-3,12,24,12,-3,-5,-2 -1,-4,-5,-3,-0,-3,-5,-4,-1 -1,-2,-4,-5,-5,-5,-4,-2,-1 0,-1,-1,-2,-2,-2,-1,-1,0"); 1367 break; 1368 } 1369 if (kernel == (KernelInfo *) NULL) 1370 return(kernel); 1371 kernel->type = type; 1372 break; 1373 } 1374 case SobelKernel: 1375 { /* Simple Sobel Kernel */ 1376 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1"); 1377 if (kernel == (KernelInfo *) NULL) 1378 return(kernel); 1379 kernel->type = type; 1380 RotateKernelInfo(kernel, args->rho); 1381 break; 1382 } 1383 case RobertsKernel: 1384 { 1385 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0"); 1386 if (kernel == (KernelInfo *) NULL) 1387 return(kernel); 1388 kernel->type = type; 1389 RotateKernelInfo(kernel, args->rho); 1390 break; 1391 } 1392 case PrewittKernel: 1393 { 1394 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1"); 1395 if (kernel == (KernelInfo *) NULL) 1396 return(kernel); 1397 kernel->type = type; 1398 RotateKernelInfo(kernel, args->rho); 1399 break; 1400 } 1401 case CompassKernel: 1402 { 1403 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1"); 1404 if (kernel == (KernelInfo *) NULL) 1405 return(kernel); 1406 kernel->type = type; 1407 RotateKernelInfo(kernel, args->rho); 1408 break; 1409 } 1410 case KirschKernel: 1411 { 1412 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3"); 1413 if (kernel == (KernelInfo *) NULL) 1414 return(kernel); 1415 kernel->type = type; 1416 RotateKernelInfo(kernel, args->rho); 1417 break; 1418 } 1419 case FreiChenKernel: 1420 /* Direction is set to be left to right positive */ 1421 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */ 1422 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */ 1423 { switch ( (int) args->rho ) { 1424 default: 1425 case 0: 1426 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1"); 1427 if (kernel == (KernelInfo *) NULL) 1428 return(kernel); 1429 kernel->type = type; 1430 kernel->values[3] = +(MagickRealType) MagickSQ2; 1431 kernel->values[5] = -(MagickRealType) MagickSQ2; 1432 CalcKernelMetaData(kernel); /* recalculate meta-data */ 1433 break; 1434 case 2: 1435 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1"); 1436 if (kernel == (KernelInfo *) NULL) 1437 return(kernel); 1438 kernel->type = type; 1439 kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2; 1440 kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2; 1441 CalcKernelMetaData(kernel); /* recalculate meta-data */ 1442 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue); 1443 break; 1444 case 10: 1445 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19"); 1446 if (kernel == (KernelInfo *) NULL) 1447 return(kernel); 1448 break; 1449 case 1: 1450 case 11: 1451 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1"); 1452 if (kernel == (KernelInfo *) NULL) 1453 return(kernel); 1454 kernel->type = type; 1455 kernel->values[3] = +(MagickRealType) MagickSQ2; 1456 kernel->values[5] = -(MagickRealType) MagickSQ2; 1457 CalcKernelMetaData(kernel); /* recalculate meta-data */ 1458 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue); 1459 break; 1460 case 12: 1461 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1"); 1462 if (kernel == (KernelInfo *) NULL) 1463 return(kernel); 1464 kernel->type = type; 1465 kernel->values[1] = +(MagickRealType) MagickSQ2; 1466 kernel->values[7] = +(MagickRealType) MagickSQ2; 1467 CalcKernelMetaData(kernel); 1468 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue); 1469 break; 1470 case 13: 1471 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2"); 1472 if (kernel == (KernelInfo *) NULL) 1473 return(kernel); 1474 kernel->type = type; 1475 kernel->values[0] = +(MagickRealType) MagickSQ2; 1476 kernel->values[8] = -(MagickRealType) MagickSQ2; 1477 CalcKernelMetaData(kernel); 1478 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue); 1479 break; 1480 case 14: 1481 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0"); 1482 if (kernel == (KernelInfo *) NULL) 1483 return(kernel); 1484 kernel->type = type; 1485 kernel->values[2] = -(MagickRealType) MagickSQ2; 1486 kernel->values[6] = +(MagickRealType) MagickSQ2; 1487 CalcKernelMetaData(kernel); 1488 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue); 1489 break; 1490 case 15: 1491 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0"); 1492 if (kernel == (KernelInfo *) NULL) 1493 return(kernel); 1494 kernel->type = type; 1495 ScaleKernelInfo(kernel, 1.0/2.0, NoValue); 1496 break; 1497 case 16: 1498 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1"); 1499 if (kernel == (KernelInfo *) NULL) 1500 return(kernel); 1501 kernel->type = type; 1502 ScaleKernelInfo(kernel, 1.0/2.0, NoValue); 1503 break; 1504 case 17: 1505 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1"); 1506 if (kernel == (KernelInfo *) NULL) 1507 return(kernel); 1508 kernel->type = type; 1509 ScaleKernelInfo(kernel, 1.0/6.0, NoValue); 1510 break; 1511 case 18: 1512 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2"); 1513 if (kernel == (KernelInfo *) NULL) 1514 return(kernel); 1515 kernel->type = type; 1516 ScaleKernelInfo(kernel, 1.0/6.0, NoValue); 1517 break; 1518 case 19: 1519 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1"); 1520 if (kernel == (KernelInfo *) NULL) 1521 return(kernel); 1522 kernel->type = type; 1523 ScaleKernelInfo(kernel, 1.0/3.0, NoValue); 1524 break; 1525 } 1526 if ( fabs(args->sigma) >= MagickEpsilon ) 1527 /* Rotate by correctly supplied 'angle' */ 1528 RotateKernelInfo(kernel, args->sigma); 1529 else if ( args->rho > 30.0 || args->rho < -30.0 ) 1530 /* Rotate by out of bounds 'type' */ 1531 RotateKernelInfo(kernel, args->rho); 1532 break; 1533 } 1534 1535 /* 1536 Boolean or Shaped Kernels 1537 */ 1538 case DiamondKernel: 1539 { 1540 if (args->rho < 1.0) 1541 kernel->width = kernel->height = 3; /* default radius = 1 */ 1542 else 1543 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 1544 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1545 1546 kernel->values=(MagickRealType *) MagickAssumeAligned( 1547 AcquireAlignedMemory(kernel->width,kernel->height* 1548 sizeof(*kernel->values))); 1549 if (kernel->values == (MagickRealType *) NULL) 1550 return(DestroyKernelInfo(kernel)); 1551 1552 /* set all kernel values within diamond area to scale given */ 1553 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 1554 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1555 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x) 1556 kernel->positive_range += kernel->values[i] = args->sigma; 1557 else 1558 kernel->values[i] = nan; 1559 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */ 1560 break; 1561 } 1562 case SquareKernel: 1563 case RectangleKernel: 1564 { double 1565 scale; 1566 if ( type == SquareKernel ) 1567 { 1568 if (args->rho < 1.0) 1569 kernel->width = kernel->height = 3; /* default radius = 1 */ 1570 else 1571 kernel->width = kernel->height = (size_t) (2*args->rho+1); 1572 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1573 scale = args->sigma; 1574 } 1575 else { 1576 /* NOTE: user defaults set in "AcquireKernelInfo()" */ 1577 if ( args->rho < 1.0 || args->sigma < 1.0 ) 1578 return(DestroyKernelInfo(kernel)); /* invalid args given */ 1579 kernel->width = (size_t)args->rho; 1580 kernel->height = (size_t)args->sigma; 1581 if ( args->xi < 0.0 || args->xi > (double)kernel->width || 1582 args->psi < 0.0 || args->psi > (double)kernel->height ) 1583 return(DestroyKernelInfo(kernel)); /* invalid args given */ 1584 kernel->x = (ssize_t) args->xi; 1585 kernel->y = (ssize_t) args->psi; 1586 scale = 1.0; 1587 } 1588 kernel->values=(MagickRealType *) MagickAssumeAligned( 1589 AcquireAlignedMemory(kernel->width,kernel->height* 1590 sizeof(*kernel->values))); 1591 if (kernel->values == (MagickRealType *) NULL) 1592 return(DestroyKernelInfo(kernel)); 1593 1594 /* set all kernel values to scale given */ 1595 u=(ssize_t) (kernel->width*kernel->height); 1596 for ( i=0; i < u; i++) 1597 kernel->values[i] = scale; 1598 kernel->minimum = kernel->maximum = scale; /* a flat shape */ 1599 kernel->positive_range = scale*u; 1600 break; 1601 } 1602 case OctagonKernel: 1603 { 1604 if (args->rho < 1.0) 1605 kernel->width = kernel->height = 5; /* default radius = 2 */ 1606 else 1607 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 1608 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1609 1610 kernel->values=(MagickRealType *) MagickAssumeAligned( 1611 AcquireAlignedMemory(kernel->width,kernel->height* 1612 sizeof(*kernel->values))); 1613 if (kernel->values == (MagickRealType *) NULL) 1614 return(DestroyKernelInfo(kernel)); 1615 1616 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 1617 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1618 if ( (labs((long) u)+labs((long) v)) <= 1619 ((long)kernel->x + (long)(kernel->x/2)) ) 1620 kernel->positive_range += kernel->values[i] = args->sigma; 1621 else 1622 kernel->values[i] = nan; 1623 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */ 1624 break; 1625 } 1626 case DiskKernel: 1627 { 1628 ssize_t 1629 limit = (ssize_t)(args->rho*args->rho); 1630 1631 if (args->rho < 0.4) /* default radius approx 4.3 */ 1632 kernel->width = kernel->height = 9L, limit = 18L; 1633 else 1634 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1; 1635 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1636 1637 kernel->values=(MagickRealType *) MagickAssumeAligned( 1638 AcquireAlignedMemory(kernel->width,kernel->height* 1639 sizeof(*kernel->values))); 1640 if (kernel->values == (MagickRealType *) NULL) 1641 return(DestroyKernelInfo(kernel)); 1642 1643 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 1644 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1645 if ((u*u+v*v) <= limit) 1646 kernel->positive_range += kernel->values[i] = args->sigma; 1647 else 1648 kernel->values[i] = nan; 1649 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */ 1650 break; 1651 } 1652 case PlusKernel: 1653 { 1654 if (args->rho < 1.0) 1655 kernel->width = kernel->height = 5; /* default radius 2 */ 1656 else 1657 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 1658 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1659 1660 kernel->values=(MagickRealType *) MagickAssumeAligned( 1661 AcquireAlignedMemory(kernel->width,kernel->height* 1662 sizeof(*kernel->values))); 1663 if (kernel->values == (MagickRealType *) NULL) 1664 return(DestroyKernelInfo(kernel)); 1665 1666 /* set all kernel values along axises to given scale */ 1667 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 1668 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1669 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan; 1670 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */ 1671 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0); 1672 break; 1673 } 1674 case CrossKernel: 1675 { 1676 if (args->rho < 1.0) 1677 kernel->width = kernel->height = 5; /* default radius 2 */ 1678 else 1679 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 1680 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1681 1682 kernel->values=(MagickRealType *) MagickAssumeAligned( 1683 AcquireAlignedMemory(kernel->width,kernel->height* 1684 sizeof(*kernel->values))); 1685 if (kernel->values == (MagickRealType *) NULL) 1686 return(DestroyKernelInfo(kernel)); 1687 1688 /* set all kernel values along axises to given scale */ 1689 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 1690 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1691 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan; 1692 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */ 1693 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0); 1694 break; 1695 } 1696 /* 1697 HitAndMiss Kernels 1698 */ 1699 case RingKernel: 1700 case PeaksKernel: 1701 { 1702 ssize_t 1703 limit1, 1704 limit2, 1705 scale; 1706 1707 if (args->rho < args->sigma) 1708 { 1709 kernel->width = ((size_t)args->sigma)*2+1; 1710 limit1 = (ssize_t)(args->rho*args->rho); 1711 limit2 = (ssize_t)(args->sigma*args->sigma); 1712 } 1713 else 1714 { 1715 kernel->width = ((size_t)args->rho)*2+1; 1716 limit1 = (ssize_t)(args->sigma*args->sigma); 1717 limit2 = (ssize_t)(args->rho*args->rho); 1718 } 1719 if ( limit2 <= 0 ) 1720 kernel->width = 7L, limit1 = 7L, limit2 = 11L; 1721 1722 kernel->height = kernel->width; 1723 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 1724 kernel->values=(MagickRealType *) MagickAssumeAligned( 1725 AcquireAlignedMemory(kernel->width,kernel->height* 1726 sizeof(*kernel->values))); 1727 if (kernel->values == (MagickRealType *) NULL) 1728 return(DestroyKernelInfo(kernel)); 1729 1730 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */ 1731 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi); 1732 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++) 1733 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 1734 { ssize_t radius=u*u+v*v; 1735 if (limit1 < radius && radius <= limit2) 1736 kernel->positive_range += kernel->values[i] = (double) scale; 1737 else 1738 kernel->values[i] = nan; 1739 } 1740 kernel->minimum = kernel->maximum = (double) scale; 1741 if ( type == PeaksKernel ) { 1742 /* set the central point in the middle */ 1743 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0; 1744 kernel->positive_range = 1.0; 1745 kernel->maximum = 1.0; 1746 } 1747 break; 1748 } 1749 case EdgesKernel: 1750 { 1751 kernel=AcquireKernelInfo("ThinSE:482"); 1752 if (kernel == (KernelInfo *) NULL) 1753 return(kernel); 1754 kernel->type = type; 1755 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */ 1756 break; 1757 } 1758 case CornersKernel: 1759 { 1760 kernel=AcquireKernelInfo("ThinSE:87"); 1761 if (kernel == (KernelInfo *) NULL) 1762 return(kernel); 1763 kernel->type = type; 1764 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */ 1765 break; 1766 } 1767 case DiagonalsKernel: 1768 { 1769 switch ( (int) args->rho ) { 1770 case 0: 1771 default: 1772 { KernelInfo 1773 *new_kernel; 1774 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-"); 1775 if (kernel == (KernelInfo *) NULL) 1776 return(kernel); 1777 kernel->type = type; 1778 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-"); 1779 if (new_kernel == (KernelInfo *) NULL) 1780 return(DestroyKernelInfo(kernel)); 1781 new_kernel->type = type; 1782 LastKernelInfo(kernel)->next = new_kernel; 1783 ExpandMirrorKernelInfo(kernel); 1784 return(kernel); 1785 } 1786 case 1: 1787 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-"); 1788 break; 1789 case 2: 1790 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-"); 1791 break; 1792 } 1793 if (kernel == (KernelInfo *) NULL) 1794 return(kernel); 1795 kernel->type = type; 1796 RotateKernelInfo(kernel, args->sigma); 1797 break; 1798 } 1799 case LineEndsKernel: 1800 { /* Kernels for finding the end of thin lines */ 1801 switch ( (int) args->rho ) { 1802 case 0: 1803 default: 1804 /* set of kernels to find all end of lines */ 1805 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>")); 1806 case 1: 1807 /* kernel for 4-connected line ends - no rotation */ 1808 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-"); 1809 break; 1810 case 2: 1811 /* kernel to add for 8-connected lines - no rotation */ 1812 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1"); 1813 break; 1814 case 3: 1815 /* kernel to add for orthogonal line ends - does not find corners */ 1816 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0"); 1817 break; 1818 case 4: 1819 /* traditional line end - fails on last T end */ 1820 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-"); 1821 break; 1822 } 1823 if (kernel == (KernelInfo *) NULL) 1824 return(kernel); 1825 kernel->type = type; 1826 RotateKernelInfo(kernel, args->sigma); 1827 break; 1828 } 1829 case LineJunctionsKernel: 1830 { /* kernels for finding the junctions of multiple lines */ 1831 switch ( (int) args->rho ) { 1832 case 0: 1833 default: 1834 /* set of kernels to find all line junctions */ 1835 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>")); 1836 case 1: 1837 /* Y Junction */ 1838 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-"); 1839 break; 1840 case 2: 1841 /* Diagonal T Junctions */ 1842 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1"); 1843 break; 1844 case 3: 1845 /* Orthogonal T Junctions */ 1846 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-"); 1847 break; 1848 case 4: 1849 /* Diagonal X Junctions */ 1850 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1"); 1851 break; 1852 case 5: 1853 /* Orthogonal X Junctions - minimal diamond kernel */ 1854 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-"); 1855 break; 1856 } 1857 if (kernel == (KernelInfo *) NULL) 1858 return(kernel); 1859 kernel->type = type; 1860 RotateKernelInfo(kernel, args->sigma); 1861 break; 1862 } 1863 case RidgesKernel: 1864 { /* Ridges - Ridge finding kernels */ 1865 KernelInfo 1866 *new_kernel; 1867 switch ( (int) args->rho ) { 1868 case 1: 1869 default: 1870 kernel=ParseKernelArray("3x1:0,1,0"); 1871 if (kernel == (KernelInfo *) NULL) 1872 return(kernel); 1873 kernel->type = type; 1874 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */ 1875 break; 1876 case 2: 1877 kernel=ParseKernelArray("4x1:0,1,1,0"); 1878 if (kernel == (KernelInfo *) NULL) 1879 return(kernel); 1880 kernel->type = type; 1881 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */ 1882 1883 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */ 1884 /* Unfortunatally we can not yet rotate a non-square kernel */ 1885 /* But then we can't flip a non-symetrical kernel either */ 1886 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0"); 1887 if (new_kernel == (KernelInfo *) NULL) 1888 return(DestroyKernelInfo(kernel)); 1889 new_kernel->type = type; 1890 LastKernelInfo(kernel)->next = new_kernel; 1891 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0"); 1892 if (new_kernel == (KernelInfo *) NULL) 1893 return(DestroyKernelInfo(kernel)); 1894 new_kernel->type = type; 1895 LastKernelInfo(kernel)->next = new_kernel; 1896 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-"); 1897 if (new_kernel == (KernelInfo *) NULL) 1898 return(DestroyKernelInfo(kernel)); 1899 new_kernel->type = type; 1900 LastKernelInfo(kernel)->next = new_kernel; 1901 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-"); 1902 if (new_kernel == (KernelInfo *) NULL) 1903 return(DestroyKernelInfo(kernel)); 1904 new_kernel->type = type; 1905 LastKernelInfo(kernel)->next = new_kernel; 1906 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0"); 1907 if (new_kernel == (KernelInfo *) NULL) 1908 return(DestroyKernelInfo(kernel)); 1909 new_kernel->type = type; 1910 LastKernelInfo(kernel)->next = new_kernel; 1911 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0"); 1912 if (new_kernel == (KernelInfo *) NULL) 1913 return(DestroyKernelInfo(kernel)); 1914 new_kernel->type = type; 1915 LastKernelInfo(kernel)->next = new_kernel; 1916 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-"); 1917 if (new_kernel == (KernelInfo *) NULL) 1918 return(DestroyKernelInfo(kernel)); 1919 new_kernel->type = type; 1920 LastKernelInfo(kernel)->next = new_kernel; 1921 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-"); 1922 if (new_kernel == (KernelInfo *) NULL) 1923 return(DestroyKernelInfo(kernel)); 1924 new_kernel->type = type; 1925 LastKernelInfo(kernel)->next = new_kernel; 1926 break; 1927 } 1928 break; 1929 } 1930 case ConvexHullKernel: 1931 { 1932 KernelInfo 1933 *new_kernel; 1934 /* first set of 8 kernels */ 1935 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0"); 1936 if (kernel == (KernelInfo *) NULL) 1937 return(kernel); 1938 kernel->type = type; 1939 ExpandRotateKernelInfo(kernel, 90.0); 1940 /* append the mirror versions too - no flip function yet */ 1941 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0"); 1942 if (new_kernel == (KernelInfo *) NULL) 1943 return(DestroyKernelInfo(kernel)); 1944 new_kernel->type = type; 1945 ExpandRotateKernelInfo(new_kernel, 90.0); 1946 LastKernelInfo(kernel)->next = new_kernel; 1947 break; 1948 } 1949 case SkeletonKernel: 1950 { 1951 switch ( (int) args->rho ) { 1952 case 1: 1953 default: 1954 /* Traditional Skeleton... 1955 ** A cyclically rotated single kernel 1956 */ 1957 kernel=AcquireKernelInfo("ThinSE:482"); 1958 if (kernel == (KernelInfo *) NULL) 1959 return(kernel); 1960 kernel->type = type; 1961 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */ 1962 break; 1963 case 2: 1964 /* HIPR Variation of the cyclic skeleton 1965 ** Corners of the traditional method made more forgiving, 1966 ** but the retain the same cyclic order. 1967 */ 1968 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;"); 1969 if (kernel == (KernelInfo *) NULL) 1970 return(kernel); 1971 if (kernel->next == (KernelInfo *) NULL) 1972 return(DestroyKernelInfo(kernel)); 1973 kernel->type = type; 1974 kernel->next->type = type; 1975 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */ 1976 break; 1977 case 3: 1978 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's 1979 ** "Connectivity-Preserving Morphological Image Thransformations" 1980 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers, 1981 ** http://www.leptonica.com/papers/conn.pdf 1982 */ 1983 kernel=AcquireKernelInfo( 1984 "ThinSE:41; ThinSE:42; ThinSE:43"); 1985 if (kernel == (KernelInfo *) NULL) 1986 return(kernel); 1987 kernel->type = type; 1988 kernel->next->type = type; 1989 kernel->next->next->type = type; 1990 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */ 1991 break; 1992 } 1993 break; 1994 } 1995 case ThinSEKernel: 1996 { /* Special kernels for general thinning, while preserving connections 1997 ** "Connectivity-Preserving Morphological Image Thransformations" 1998 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers, 1999 ** http://www.leptonica.com/papers/conn.pdf 2000 ** And 2001 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html 2002 ** 2003 ** Note kernels do not specify the origin pixel, allowing them 2004 ** to be used for both thickening and thinning operations. 2005 */ 2006 switch ( (int) args->rho ) { 2007 /* SE for 4-connected thinning */ 2008 case 41: /* SE_4_1 */ 2009 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1"); 2010 break; 2011 case 42: /* SE_4_2 */ 2012 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-"); 2013 break; 2014 case 43: /* SE_4_3 */ 2015 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1"); 2016 break; 2017 case 44: /* SE_4_4 */ 2018 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-"); 2019 break; 2020 case 45: /* SE_4_5 */ 2021 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-"); 2022 break; 2023 case 46: /* SE_4_6 */ 2024 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1"); 2025 break; 2026 case 47: /* SE_4_7 */ 2027 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-"); 2028 break; 2029 case 48: /* SE_4_8 */ 2030 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1"); 2031 break; 2032 case 49: /* SE_4_9 */ 2033 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1"); 2034 break; 2035 /* SE for 8-connected thinning - negatives of the above */ 2036 case 81: /* SE_8_0 */ 2037 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-"); 2038 break; 2039 case 82: /* SE_8_2 */ 2040 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-"); 2041 break; 2042 case 83: /* SE_8_3 */ 2043 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-"); 2044 break; 2045 case 84: /* SE_8_4 */ 2046 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-"); 2047 break; 2048 case 85: /* SE_8_5 */ 2049 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-"); 2050 break; 2051 case 86: /* SE_8_6 */ 2052 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1"); 2053 break; 2054 case 87: /* SE_8_7 */ 2055 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-"); 2056 break; 2057 case 88: /* SE_8_8 */ 2058 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-"); 2059 break; 2060 case 89: /* SE_8_9 */ 2061 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-"); 2062 break; 2063 /* Special combined SE kernels */ 2064 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */ 2065 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-"); 2066 break; 2067 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */ 2068 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-"); 2069 break; 2070 case 481: /* SE_48_1 - General Connected Corner Kernel */ 2071 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-"); 2072 break; 2073 default: 2074 case 482: /* SE_48_2 - General Edge Kernel */ 2075 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1"); 2076 break; 2077 } 2078 if (kernel == (KernelInfo *) NULL) 2079 return(kernel); 2080 kernel->type = type; 2081 RotateKernelInfo(kernel, args->sigma); 2082 break; 2083 } 2084 /* 2085 Distance Measuring Kernels 2086 */ 2087 case ChebyshevKernel: 2088 { 2089 if (args->rho < 1.0) 2090 kernel->width = kernel->height = 3; /* default radius = 1 */ 2091 else 2092 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 2093 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 2094 2095 kernel->values=(MagickRealType *) MagickAssumeAligned( 2096 AcquireAlignedMemory(kernel->width,kernel->height* 2097 sizeof(*kernel->values))); 2098 if (kernel->values == (MagickRealType *) NULL) 2099 return(DestroyKernelInfo(kernel)); 2100 2101 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 2102 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 2103 kernel->positive_range += ( kernel->values[i] = 2104 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) ); 2105 kernel->maximum = kernel->values[0]; 2106 break; 2107 } 2108 case ManhattanKernel: 2109 { 2110 if (args->rho < 1.0) 2111 kernel->width = kernel->height = 3; /* default radius = 1 */ 2112 else 2113 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 2114 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 2115 2116 kernel->values=(MagickRealType *) MagickAssumeAligned( 2117 AcquireAlignedMemory(kernel->width,kernel->height* 2118 sizeof(*kernel->values))); 2119 if (kernel->values == (MagickRealType *) NULL) 2120 return(DestroyKernelInfo(kernel)); 2121 2122 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 2123 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 2124 kernel->positive_range += ( kernel->values[i] = 2125 args->sigma*(labs((long) u)+labs((long) v)) ); 2126 kernel->maximum = kernel->values[0]; 2127 break; 2128 } 2129 case OctagonalKernel: 2130 { 2131 if (args->rho < 2.0) 2132 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */ 2133 else 2134 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 2135 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 2136 2137 kernel->values=(MagickRealType *) MagickAssumeAligned( 2138 AcquireAlignedMemory(kernel->width,kernel->height* 2139 sizeof(*kernel->values))); 2140 if (kernel->values == (MagickRealType *) NULL) 2141 return(DestroyKernelInfo(kernel)); 2142 2143 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 2144 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 2145 { 2146 double 2147 r1 = MagickMax(fabs((double)u),fabs((double)v)), 2148 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5); 2149 kernel->positive_range += kernel->values[i] = 2150 args->sigma*MagickMax(r1,r2); 2151 } 2152 kernel->maximum = kernel->values[0]; 2153 break; 2154 } 2155 case EuclideanKernel: 2156 { 2157 if (args->rho < 1.0) 2158 kernel->width = kernel->height = 3; /* default radius = 1 */ 2159 else 2160 kernel->width = kernel->height = ((size_t)args->rho)*2+1; 2161 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2; 2162 2163 kernel->values=(MagickRealType *) MagickAssumeAligned( 2164 AcquireAlignedMemory(kernel->width,kernel->height* 2165 sizeof(*kernel->values))); 2166 if (kernel->values == (MagickRealType *) NULL) 2167 return(DestroyKernelInfo(kernel)); 2168 2169 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++) 2170 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++) 2171 kernel->positive_range += ( kernel->values[i] = 2172 args->sigma*sqrt((double)(u*u+v*v)) ); 2173 kernel->maximum = kernel->values[0]; 2174 break; 2175 } 2176 default: 2177 { 2178 /* No-Op Kernel - Basically just a single pixel on its own */ 2179 kernel=ParseKernelArray("1:1"); 2180 if (kernel == (KernelInfo *) NULL) 2181 return(kernel); 2182 kernel->type = UndefinedKernel; 2183 break; 2184 } 2185 break; 2186 } 2187 return(kernel); 2188} 2189 2190/* 2191%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2192% % 2193% % 2194% % 2195% C l o n e K e r n e l I n f o % 2196% % 2197% % 2198% % 2199%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2200% 2201% CloneKernelInfo() creates a new clone of the given Kernel List so that its 2202% can be modified without effecting the original. The cloned kernel should 2203% be destroyed using DestoryKernelInfo() when no longer needed. 2204% 2205% The format of the CloneKernelInfo method is: 2206% 2207% KernelInfo *CloneKernelInfo(const KernelInfo *kernel) 2208% 2209% A description of each parameter follows: 2210% 2211% o kernel: the Morphology/Convolution kernel to be cloned 2212% 2213*/ 2214MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel) 2215{ 2216 register ssize_t 2217 i; 2218 2219 KernelInfo 2220 *new_kernel; 2221 2222 assert(kernel != (KernelInfo *) NULL); 2223 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel)); 2224 if (new_kernel == (KernelInfo *) NULL) 2225 return(new_kernel); 2226 *new_kernel=(*kernel); /* copy values in structure */ 2227 2228 /* replace the values with a copy of the values */ 2229 new_kernel->values=(MagickRealType *) MagickAssumeAligned( 2230 AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values))); 2231 if (new_kernel->values == (MagickRealType *) NULL) 2232 return(DestroyKernelInfo(new_kernel)); 2233 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++) 2234 new_kernel->values[i]=kernel->values[i]; 2235 2236 /* Also clone the next kernel in the kernel list */ 2237 if ( kernel->next != (KernelInfo *) NULL ) { 2238 new_kernel->next = CloneKernelInfo(kernel->next); 2239 if ( new_kernel->next == (KernelInfo *) NULL ) 2240 return(DestroyKernelInfo(new_kernel)); 2241 } 2242 2243 return(new_kernel); 2244} 2245 2246/* 2247%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2248% % 2249% % 2250% % 2251% D e s t r o y K e r n e l I n f o % 2252% % 2253% % 2254% % 2255%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2256% 2257% DestroyKernelInfo() frees the memory used by a Convolution/Morphology 2258% kernel. 2259% 2260% The format of the DestroyKernelInfo method is: 2261% 2262% KernelInfo *DestroyKernelInfo(KernelInfo *kernel) 2263% 2264% A description of each parameter follows: 2265% 2266% o kernel: the Morphology/Convolution kernel to be destroyed 2267% 2268*/ 2269MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel) 2270{ 2271 assert(kernel != (KernelInfo *) NULL); 2272 if ( kernel->next != (KernelInfo *) NULL ) 2273 kernel->next=DestroyKernelInfo(kernel->next); 2274 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values); 2275 kernel=(KernelInfo *) RelinquishMagickMemory(kernel); 2276 return(kernel); 2277} 2278 2279/* 2280%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2281% % 2282% % 2283% % 2284+ E x p a n d M i r r o r K e r n e l I n f o % 2285% % 2286% % 2287% % 2288%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2289% 2290% ExpandMirrorKernelInfo() takes a single kernel, and expands it into a 2291% sequence of 90-degree rotated kernels but providing a reflected 180 2292% rotatation, before the -/+ 90-degree rotations. 2293% 2294% This special rotation order produces a better, more symetrical thinning of 2295% objects. 2296% 2297% The format of the ExpandMirrorKernelInfo method is: 2298% 2299% void ExpandMirrorKernelInfo(KernelInfo *kernel) 2300% 2301% A description of each parameter follows: 2302% 2303% o kernel: the Morphology/Convolution kernel 2304% 2305% This function is only internel to this module, as it is not finalized, 2306% especially with regard to non-orthogonal angles, and rotation of larger 2307% 2D kernels. 2308*/ 2309 2310#if 0 2311static void FlopKernelInfo(KernelInfo *kernel) 2312 { /* Do a Flop by reversing each row. */ 2313 size_t 2314 y; 2315 register ssize_t 2316 x,r; 2317 register double 2318 *k,t; 2319 2320 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width) 2321 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--) 2322 t=k[x], k[x]=k[r], k[r]=t; 2323 2324 kernel->x = kernel->width - kernel->x - 1; 2325 angle = fmod(angle+180.0, 360.0); 2326 } 2327#endif 2328 2329static void ExpandMirrorKernelInfo(KernelInfo *kernel) 2330{ 2331 KernelInfo 2332 *clone, 2333 *last; 2334 2335 last = kernel; 2336 2337 clone = CloneKernelInfo(last); 2338 RotateKernelInfo(clone, 180); /* flip */ 2339 LastKernelInfo(last)->next = clone; 2340 last = clone; 2341 2342 clone = CloneKernelInfo(last); 2343 RotateKernelInfo(clone, 90); /* transpose */ 2344 LastKernelInfo(last)->next = clone; 2345 last = clone; 2346 2347 clone = CloneKernelInfo(last); 2348 RotateKernelInfo(clone, 180); /* flop */ 2349 LastKernelInfo(last)->next = clone; 2350 2351 return; 2352} 2353 2354/* 2355%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2356% % 2357% % 2358% % 2359+ E x p a n d R o t a t e K e r n e l I n f o % 2360% % 2361% % 2362% % 2363%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2364% 2365% ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating 2366% incrementally by the angle given, until the kernel repeats. 2367% 2368% WARNING: 45 degree rotations only works for 3x3 kernels. 2369% While 90 degree roatations only works for linear and square kernels 2370% 2371% The format of the ExpandRotateKernelInfo method is: 2372% 2373% void ExpandRotateKernelInfo(KernelInfo *kernel, double angle) 2374% 2375% A description of each parameter follows: 2376% 2377% o kernel: the Morphology/Convolution kernel 2378% 2379% o angle: angle to rotate in degrees 2380% 2381% This function is only internel to this module, as it is not finalized, 2382% especially with regard to non-orthogonal angles, and rotation of larger 2383% 2D kernels. 2384*/ 2385 2386/* Internal Routine - Return true if two kernels are the same */ 2387static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1, 2388 const KernelInfo *kernel2) 2389{ 2390 register size_t 2391 i; 2392 2393 /* check size and origin location */ 2394 if ( kernel1->width != kernel2->width 2395 || kernel1->height != kernel2->height 2396 || kernel1->x != kernel2->x 2397 || kernel1->y != kernel2->y ) 2398 return MagickFalse; 2399 2400 /* check actual kernel values */ 2401 for (i=0; i < (kernel1->width*kernel1->height); i++) { 2402 /* Test for Nan equivalence */ 2403 if ( IfNaN(kernel1->values[i]) && !IfNaN(kernel2->values[i]) ) 2404 return MagickFalse; 2405 if ( IfNaN(kernel2->values[i]) && !IfNaN(kernel1->values[i]) ) 2406 return MagickFalse; 2407 /* Test actual values are equivalent */ 2408 if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon ) 2409 return MagickFalse; 2410 } 2411 2412 return MagickTrue; 2413} 2414 2415static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle) 2416{ 2417 KernelInfo 2418 *clone, 2419 *last; 2420 2421 last = kernel; 2422 while(1) { 2423 clone = CloneKernelInfo(last); 2424 RotateKernelInfo(clone, angle); 2425 if ( SameKernelInfo(kernel, clone) == MagickTrue ) 2426 break; 2427 LastKernelInfo(last)->next = clone; 2428 last = clone; 2429 } 2430 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */ 2431 return; 2432} 2433 2434/* 2435%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2436% % 2437% % 2438% % 2439+ C a l c M e t a K e r n a l I n f o % 2440% % 2441% % 2442% % 2443%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2444% 2445% CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only, 2446% using the kernel values. This should only ne used if it is not possible to 2447% calculate that meta-data in some easier way. 2448% 2449% It is important that the meta-data is correct before ScaleKernelInfo() is 2450% used to perform kernel normalization. 2451% 2452% The format of the CalcKernelMetaData method is: 2453% 2454% void CalcKernelMetaData(KernelInfo *kernel, const double scale ) 2455% 2456% A description of each parameter follows: 2457% 2458% o kernel: the Morphology/Convolution kernel to modify 2459% 2460% WARNING: Minimum and Maximum values are assumed to include zero, even if 2461% zero is not part of the kernel (as in Gaussian Derived kernels). This 2462% however is not true for flat-shaped morphological kernels. 2463% 2464% WARNING: Only the specific kernel pointed to is modified, not a list of 2465% multiple kernels. 2466% 2467% This is an internal function and not expected to be useful outside this 2468% module. This could change however. 2469*/ 2470static void CalcKernelMetaData(KernelInfo *kernel) 2471{ 2472 register size_t 2473 i; 2474 2475 kernel->minimum = kernel->maximum = 0.0; 2476 kernel->negative_range = kernel->positive_range = 0.0; 2477 for (i=0; i < (kernel->width*kernel->height); i++) 2478 { 2479 if ( fabs(kernel->values[i]) < MagickEpsilon ) 2480 kernel->values[i] = 0.0; 2481 ( kernel->values[i] < 0) 2482 ? ( kernel->negative_range += kernel->values[i] ) 2483 : ( kernel->positive_range += kernel->values[i] ); 2484 Minimize(kernel->minimum, kernel->values[i]); 2485 Maximize(kernel->maximum, kernel->values[i]); 2486 } 2487 2488 return; 2489} 2490 2491/* 2492%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2493% % 2494% % 2495% % 2496% M o r p h o l o g y A p p l y % 2497% % 2498% % 2499% % 2500%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2501% 2502% MorphologyApply() applies a morphological method, multiple times using 2503% a list of multiple kernels. This is the method that should be called by 2504% other 'operators' that internally use morphology operations as part of 2505% their processing. 2506% 2507% It is basically equivalent to as MorphologyImage() (see below) but without 2508% any user controls. This allows internel programs to use this method to 2509% perform a specific task without possible interference by any API user 2510% supplied settings. 2511% 2512% It is MorphologyImage() task to extract any such user controls, and 2513% pass them to this function for processing. 2514% 2515% More specifically all given kernels should already be scaled, normalised, 2516% and blended appropriatally before being parred to this routine. The 2517% appropriate bias, and compose (typically 'UndefinedComposeOp') given. 2518% 2519% The format of the MorphologyApply method is: 2520% 2521% Image *MorphologyApply(const Image *image,MorphologyMethod method, 2522% const ssize_t iterations,const KernelInfo *kernel, 2523% const CompositeMethod compose,const double bias, 2524% ExceptionInfo *exception) 2525% 2526% A description of each parameter follows: 2527% 2528% o image: the source image 2529% 2530% o method: the morphology method to be applied. 2531% 2532% o iterations: apply the operation this many times (or no change). 2533% A value of -1 means loop until no change found. 2534% How this is applied may depend on the morphology method. 2535% Typically this is a value of 1. 2536% 2537% o channel: the channel type. 2538% 2539% o kernel: An array of double representing the morphology kernel. 2540% 2541% o compose: How to handle or merge multi-kernel results. 2542% If 'UndefinedCompositeOp' use default for the Morphology method. 2543% If 'NoCompositeOp' force image to be re-iterated by each kernel. 2544% Otherwise merge the results using the compose method given. 2545% 2546% o bias: Convolution Output Bias. 2547% 2548% o exception: return any errors or warnings in this structure. 2549% 2550*/ 2551static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image, 2552 const MorphologyMethod method,const KernelInfo *kernel,const double bias, 2553 ExceptionInfo *exception) 2554{ 2555#define MorphologyTag "Morphology/Image" 2556 2557 CacheView 2558 *image_view, 2559 *morphology_view; 2560 2561 OffsetInfo 2562 offset; 2563 2564 register ssize_t 2565 i; 2566 2567 ssize_t 2568 y; 2569 2570 size_t 2571 *changes, 2572 changed, 2573 width; 2574 2575 MagickBooleanType 2576 status; 2577 2578 MagickOffsetType 2579 progress; 2580 2581 assert(image != (Image *) NULL); 2582 assert(image->signature == MagickSignature); 2583 assert(morphology_image != (Image *) NULL); 2584 assert(morphology_image->signature == MagickSignature); 2585 assert(kernel != (KernelInfo *) NULL); 2586 assert(kernel->signature == MagickSignature); 2587 assert(exception != (ExceptionInfo *) NULL); 2588 assert(exception->signature == MagickSignature); 2589 status=MagickTrue; 2590 progress=0; 2591 image_view=AcquireVirtualCacheView(image,exception); 2592 morphology_view=AcquireAuthenticCacheView(morphology_image,exception); 2593 width=image->columns+kernel->width-1; 2594 switch (method) 2595 { 2596 case ConvolveMorphology: 2597 case DilateMorphology: 2598 case DilateIntensityMorphology: 2599 case IterativeDistanceMorphology: 2600 { 2601 /* 2602 Kernel needs to used with reflection about origin. 2603 */ 2604 offset.x=(ssize_t) kernel->width-kernel->x-1; 2605 offset.y=(ssize_t) kernel->height-kernel->y-1; 2606 break; 2607 } 2608 case ErodeMorphology: 2609 case ErodeIntensityMorphology: 2610 case HitAndMissMorphology: 2611 case ThinningMorphology: 2612 case ThickenMorphology: 2613 { 2614 offset.x=kernel->x; 2615 offset.y=kernel->y; 2616 break; 2617 } 2618 default: 2619 { 2620 assert("Not a Primitive Morphology Method" != (char *) NULL); 2621 break; 2622 } 2623 } 2624 changed=0; 2625 changes=(size_t *) AcquireQuantumMemory(GetOpenMPMaximumThreads(), 2626 sizeof(*changes)); 2627 if (changes == (size_t *) NULL) 2628 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed"); 2629 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++) 2630 changes[i]=0; 2631 if ((method == ConvolveMorphology) && (kernel->width == 1)) 2632 { 2633 const int 2634 id = GetOpenMPThreadId(); 2635 2636 register ssize_t 2637 x; 2638 2639 /* 2640 Special handling (for speed) of vertical (blur) kernels. This performs 2641 its handling in columns rather than in rows. This is only done 2642 for convolve as it is the only method that generates very large 1-D 2643 vertical kernels (such as a 'BlurKernel') 2644 */ 2645#if defined(MAGICKCORE_OPENMP_SUPPORT) 2646 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 2647 magick_threads(image,morphology_image,image->columns,1) 2648#endif 2649 for (x=0; x < (ssize_t) image->columns; x++) 2650 { 2651 register const Quantum 2652 *restrict p; 2653 2654 register Quantum 2655 *restrict q; 2656 2657 register ssize_t 2658 y; 2659 2660 ssize_t 2661 center; 2662 2663 if (status == MagickFalse) 2664 continue; 2665 p=GetCacheViewVirtualPixels(image_view,x,-offset.y,1,image->rows+ 2666 kernel->height-1,exception); 2667 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1, 2668 morphology_image->rows,exception); 2669 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL)) 2670 { 2671 status=MagickFalse; 2672 continue; 2673 } 2674 center=(ssize_t) GetPixelChannels(image)*offset.y; 2675 for (y=0; y < (ssize_t) image->rows; y++) 2676 { 2677 register ssize_t 2678 i; 2679 2680 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 2681 { 2682 double 2683 alpha, 2684 pixel; 2685 2686 PixelChannel 2687 channel; 2688 2689 PixelTrait 2690 morphology_traits, 2691 traits; 2692 2693 register const MagickRealType 2694 *restrict k; 2695 2696 register const Quantum 2697 *restrict pixels; 2698 2699 register ssize_t 2700 u; 2701 2702 ssize_t 2703 v; 2704 2705 channel=GetPixelChannelChannel(image,i); 2706 traits=GetPixelChannelTraits(image,channel); 2707 morphology_traits=GetPixelChannelTraits(morphology_image,channel); 2708 if ((traits == UndefinedPixelTrait) || 2709 (morphology_traits == UndefinedPixelTrait)) 2710 continue; 2711 if (((morphology_traits & CopyPixelTrait) != 0) || 2712 (GetPixelReadMask(image,p+center) == 0)) 2713 { 2714 SetPixelChannel(morphology_image,channel,p[center+i],q); 2715 continue; 2716 } 2717 k=(&kernel->values[kernel->width*kernel->height-1]); 2718 pixels=p; 2719 pixel=bias; 2720 if ((morphology_traits & BlendPixelTrait) == 0) 2721 { 2722 /* 2723 No alpha blending. 2724 */ 2725 for (v=0; v < (ssize_t) kernel->height; v++) 2726 { 2727 for (u=0; u < (ssize_t) kernel->width; u++) 2728 { 2729 if (IfNaN(*k) == MagickFalse) 2730 pixel+=(*k)*pixels[i]; 2731 k--; 2732 pixels+=GetPixelChannels(image); 2733 } 2734 } 2735 if (fabs(pixel-p[center+i]) > MagickEpsilon) 2736 changes[id]++; 2737 SetPixelChannel(morphology_image,channel,ClampToQuantum(pixel), 2738 q); 2739 continue; 2740 } 2741 /* 2742 Alpha blending. 2743 */ 2744 for (v=0; v < (ssize_t) kernel->height; v++) 2745 { 2746 for (u=0; u < (ssize_t) kernel->width; u++) 2747 { 2748 if (IfNaN(*k) == MagickFalse) 2749 { 2750 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels)); 2751 pixel+=(*k)*alpha*pixels[i]; 2752 } 2753 k--; 2754 pixels+=GetPixelChannels(image); 2755 } 2756 } 2757 if (fabs(pixel-p[center+i]) > MagickEpsilon) 2758 changes[id]++; 2759 SetPixelChannel(morphology_image,channel,ClampToQuantum(pixel),q); 2760 } 2761 p+=GetPixelChannels(image); 2762 q+=GetPixelChannels(morphology_image); 2763 } 2764 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse) 2765 status=MagickFalse; 2766 if (image->progress_monitor != (MagickProgressMonitor) NULL) 2767 { 2768 MagickBooleanType 2769 proceed; 2770 2771#if defined(MAGICKCORE_OPENMP_SUPPORT) 2772 #pragma omp critical (MagickCore_MorphologyPrimitive) 2773#endif 2774 proceed=SetImageProgress(image,MorphologyTag,progress++, 2775 image->rows); 2776 if (proceed == MagickFalse) 2777 status=MagickFalse; 2778 } 2779 } 2780 morphology_image->type=image->type; 2781 morphology_view=DestroyCacheView(morphology_view); 2782 image_view=DestroyCacheView(image_view); 2783 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++) 2784 changed+=changes[i]; 2785 changes=(size_t *) RelinquishMagickMemory(changes); 2786 return(status ? (ssize_t) changed : 0); 2787 } 2788 /* 2789 Normal handling of horizontal or rectangular kernels (row by row). 2790 */ 2791#if defined(MAGICKCORE_OPENMP_SUPPORT) 2792 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 2793 magick_threads(image,morphology_image,image->rows,1) 2794#endif 2795 for (y=0; y < (ssize_t) image->rows; y++) 2796 { 2797 const int 2798 id = GetOpenMPThreadId(); 2799 2800 register const Quantum 2801 *restrict p; 2802 2803 register Quantum 2804 *restrict q; 2805 2806 register ssize_t 2807 x; 2808 2809 ssize_t 2810 center; 2811 2812 if (status == MagickFalse) 2813 continue; 2814 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width, 2815 kernel->height,exception); 2816 q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns, 2817 1,exception); 2818 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL)) 2819 { 2820 status=MagickFalse; 2821 continue; 2822 } 2823 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+ 2824 GetPixelChannels(image)*offset.x); 2825 for (x=0; x < (ssize_t) image->columns; x++) 2826 { 2827 register ssize_t 2828 i; 2829 2830 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 2831 { 2832 double 2833 alpha, 2834 maximum, 2835 minimum, 2836 pixel; 2837 2838 PixelChannel 2839 channel; 2840 2841 PixelTrait 2842 morphology_traits, 2843 traits; 2844 2845 register const MagickRealType 2846 *restrict k; 2847 2848 register const Quantum 2849 *restrict pixels; 2850 2851 register ssize_t 2852 u; 2853 2854 ssize_t 2855 v; 2856 2857 channel=GetPixelChannelChannel(image,i); 2858 traits=GetPixelChannelTraits(image,channel); 2859 morphology_traits=GetPixelChannelTraits(morphology_image,channel); 2860 if ((traits == UndefinedPixelTrait) || 2861 (morphology_traits == UndefinedPixelTrait)) 2862 continue; 2863 if (((morphology_traits & CopyPixelTrait) != 0) || 2864 (GetPixelReadMask(image,p+center) == 0)) 2865 { 2866 SetPixelChannel(morphology_image,channel,p[center+i],q); 2867 continue; 2868 } 2869 pixels=p; 2870 maximum=0.0; 2871 minimum=(double) QuantumRange; 2872 switch (method) 2873 { 2874 case ConvolveMorphology: pixel=bias; break; 2875 case HitAndMissMorphology: pixel=(double) QuantumRange; break; 2876 case ThinningMorphology: pixel=(double) QuantumRange; break; 2877 case ThickenMorphology: pixel=(double) QuantumRange; break; 2878 case ErodeMorphology: pixel=(double) QuantumRange; break; 2879 case DilateMorphology: pixel=0.0; break; 2880 case ErodeIntensityMorphology: 2881 case DilateIntensityMorphology: 2882 case IterativeDistanceMorphology: 2883 { 2884 pixel=(double) p[center+i]; 2885 break; 2886 } 2887 default: pixel=0; break; 2888 } 2889 switch (method) 2890 { 2891 case ConvolveMorphology: 2892 { 2893 /* 2894 Weighted Average of pixels using reflected kernel 2895 2896 For correct working of this operation for asymetrical 2897 kernels, the kernel needs to be applied in its reflected form. 2898 That is its values needs to be reversed. 2899 2900 Correlation is actually the same as this but without reflecting 2901 the kernel, and thus 'lower-level' that Convolution. However 2902 as Convolution is the more common method used, and it does not 2903 really cost us much in terms of processing to use a reflected 2904 kernel, so it is Convolution that is implemented. 2905 2906 Correlation will have its kernel reflected before calling 2907 this function to do a Convolve. 2908 2909 For more details of Correlation vs Convolution see 2910 http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf 2911 */ 2912 k=(&kernel->values[kernel->width*kernel->height-1]); 2913 if ((morphology_traits & BlendPixelTrait) == 0) 2914 { 2915 /* 2916 No alpha blending. 2917 */ 2918 for (v=0; v < (ssize_t) kernel->height; v++) 2919 { 2920 for (u=0; u < (ssize_t) kernel->width; u++) 2921 { 2922 if (IfNaN(*k) == MagickFalse) 2923 pixel+=(*k)*pixels[i]; 2924 k--; 2925 pixels+=GetPixelChannels(image); 2926 } 2927 pixels+=(image->columns-1)*GetPixelChannels(image); 2928 } 2929 break; 2930 } 2931 /* 2932 Alpha blending. 2933 */ 2934 for (v=0; v < (ssize_t) kernel->height; v++) 2935 { 2936 for (u=0; u < (ssize_t) kernel->width; u++) 2937 { 2938 if (IfNaN(*k) == MagickFalse) 2939 { 2940 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels)); 2941 pixel+=(*k)*alpha*pixels[i]; 2942 } 2943 k--; 2944 pixels+=GetPixelChannels(image); 2945 } 2946 pixels+=(image->columns-1)*GetPixelChannels(image); 2947 } 2948 break; 2949 } 2950 case ErodeMorphology: 2951 { 2952 /* 2953 Minimum value within kernel neighbourhood. 2954 2955 The kernel is not reflected for this operation. In normal 2956 Greyscale Morphology, the kernel value should be added 2957 to the real value, this is currently not done, due to the 2958 nature of the boolean kernels being used. 2959 */ 2960 k=kernel->values; 2961 for (v=0; v < (ssize_t) kernel->height; v++) 2962 { 2963 for (u=0; u < (ssize_t) kernel->width; u++) 2964 { 2965 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5)) 2966 { 2967 if ((double) pixels[i] < pixel) 2968 pixel=(double) pixels[i]; 2969 } 2970 k++; 2971 pixels+=GetPixelChannels(image); 2972 } 2973 pixels+=(image->columns-1)*GetPixelChannels(image); 2974 } 2975 break; 2976 } 2977 case DilateMorphology: 2978 { 2979 /* 2980 Maximum value within kernel neighbourhood. 2981 2982 For correct working of this operation for asymetrical kernels, 2983 the kernel needs to be applied in its reflected form. That is 2984 its values needs to be reversed. 2985 2986 In normal Greyscale Morphology, the kernel value should be 2987 added to the real value, this is currently not done, due to the 2988 nature of the boolean kernels being used. 2989 */ 2990 k=(&kernel->values[kernel->width*kernel->height-1]); 2991 for (v=0; v < (ssize_t) kernel->height; v++) 2992 { 2993 for (u=0; u < (ssize_t) kernel->width; u++) 2994 { 2995 if ((IfNaN(*k) == MagickFalse) && (*k > 0.5)) 2996 { 2997 if ((double) pixels[i] > pixel) 2998 pixel=(double) pixels[i]; 2999 } 3000 k--; 3001 pixels+=GetPixelChannels(image); 3002 } 3003 pixels+=(image->columns-1)*GetPixelChannels(image); 3004 } 3005 break; 3006 } 3007 case HitAndMissMorphology: 3008 case ThinningMorphology: 3009 case ThickenMorphology: 3010 { 3011 /* 3012 Minimum of foreground pixel minus maxumum of background pixels. 3013 3014 The kernel is not reflected for this operation, and consists 3015 of both foreground and background pixel neighbourhoods, 0.0 for 3016 background, and 1.0 for foreground with either Nan or 0.5 values 3017 for don't care. 3018 3019 This never produces a meaningless negative result. Such results 3020 cause Thinning/Thicken to not work correctly when used against a 3021 greyscale image. 3022 */ 3023 k=kernel->values; 3024 for (v=0; v < (ssize_t) kernel->height; v++) 3025 { 3026 for (u=0; u < (ssize_t) kernel->width; u++) 3027 { 3028 if (IfNaN(*k) == MagickFalse) 3029 { 3030 if (*k > 0.7) 3031 { 3032 if ((double) pixels[i] < pixel) 3033 pixel=(double) pixels[i]; 3034 } 3035 else 3036 if (*k < 0.3) 3037 { 3038 if ((double) pixels[i] > maximum) 3039 maximum=(double) pixels[i]; 3040 } 3041 } 3042 k++; 3043 pixels+=GetPixelChannels(image); 3044 } 3045 pixels+=(image->columns-1)*GetPixelChannels(image); 3046 } 3047 pixel-=maximum; 3048 if (pixel < 0.0) 3049 pixel=0.0; 3050 if (method == ThinningMorphology) 3051 pixel=(double) p[center+i]-pixel; 3052 else 3053 if (method == ThickenMorphology) 3054 pixel+=(double) p[center+i]+pixel; 3055 break; 3056 } 3057 case ErodeIntensityMorphology: 3058 { 3059 /* 3060 Select pixel with minimum intensity within kernel neighbourhood. 3061 3062 The kernel is not reflected for this operation. 3063 */ 3064 k=kernel->values; 3065 for (v=0; v < (ssize_t) kernel->height; v++) 3066 { 3067 for (u=0; u < (ssize_t) kernel->width; u++) 3068 { 3069 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5)) 3070 { 3071 if (GetPixelIntensity(image,pixels) < minimum) 3072 { 3073 pixel=(double) pixels[i]; 3074 minimum=GetPixelIntensity(image,pixels); 3075 } 3076 } 3077 k++; 3078 pixels+=GetPixelChannels(image); 3079 } 3080 pixels+=(image->columns-1)*GetPixelChannels(image); 3081 } 3082 break; 3083 } 3084 case DilateIntensityMorphology: 3085 { 3086 /* 3087 Select pixel with maximum intensity within kernel neighbourhood. 3088 3089 The kernel is not reflected for this operation. 3090 */ 3091 k=(&kernel->values[kernel->width*kernel->height-1]); 3092 for (v=0; v < (ssize_t) kernel->height; v++) 3093 { 3094 for (u=0; u < (ssize_t) kernel->width; u++) 3095 { 3096 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5)) 3097 { 3098 if (GetPixelIntensity(image,pixels) > maximum) 3099 { 3100 pixel=(double) pixels[i]; 3101 maximum=GetPixelIntensity(image,pixels); 3102 } 3103 } 3104 k--; 3105 pixels+=GetPixelChannels(image); 3106 } 3107 pixels+=(image->columns-1)*GetPixelChannels(image); 3108 } 3109 break; 3110 } 3111 case IterativeDistanceMorphology: 3112 { 3113 /* 3114 Compute th iterative distance from black edge of a white image 3115 shape. Essentually white values are decreased to the smallest 3116 'distance from edge' it can find. 3117 3118 It works by adding kernel values to the neighbourhood, and and 3119 select the minimum value found. The kernel is rotated before 3120 use, so kernel distances match resulting distances, when a user 3121 provided asymmetric kernel is applied. 3122 3123 This code is nearly identical to True GrayScale Morphology but 3124 not quite. 3125 3126 GreyDilate Kernel values added, maximum value found Kernel is 3127 rotated before use. 3128 3129 GrayErode: Kernel values subtracted and minimum value found No 3130 kernel rotation used. 3131 3132 Note the the Iterative Distance method is essentially a 3133 GrayErode, but with negative kernel values, and kernel rotation 3134 applied. 3135 */ 3136 k=(&kernel->values[kernel->width*kernel->height-1]); 3137 for (v=0; v < (ssize_t) kernel->height; v++) 3138 { 3139 for (u=0; u < (ssize_t) kernel->width; u++) 3140 { 3141 if (IfNaN(*k) == MagickFalse) 3142 { 3143 if ((pixels[i]+(*k)) < pixel) 3144 pixel=(double) pixels[i]+(*k); 3145 } 3146 k--; 3147 pixels+=GetPixelChannels(image); 3148 } 3149 pixels+=(image->columns-1)*GetPixelChannels(image); 3150 } 3151 break; 3152 } 3153 case UndefinedMorphology: 3154 default: 3155 break; 3156 } 3157 if (fabs(pixel-p[center+i]) > MagickEpsilon) 3158 changes[id]++; 3159 SetPixelChannel(morphology_image,channel,ClampToQuantum(pixel),q); 3160 } 3161 p+=GetPixelChannels(image); 3162 q+=GetPixelChannels(morphology_image); 3163 } 3164 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse) 3165 status=MagickFalse; 3166 if (image->progress_monitor != (MagickProgressMonitor) NULL) 3167 { 3168 MagickBooleanType 3169 proceed; 3170 3171#if defined(MAGICKCORE_OPENMP_SUPPORT) 3172 #pragma omp critical (MagickCore_MorphologyPrimitive) 3173#endif 3174 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows); 3175 if (proceed == MagickFalse) 3176 status=MagickFalse; 3177 } 3178 } 3179 morphology_view=DestroyCacheView(morphology_view); 3180 image_view=DestroyCacheView(image_view); 3181 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++) 3182 changed+=changes[i]; 3183 changes=(size_t *) RelinquishMagickMemory(changes); 3184 return(status ? (ssize_t) changed : -1); 3185} 3186 3187/* 3188 This is almost identical to the MorphologyPrimative() function above, but 3189 applies the primitive directly to the actual image using two passes, once in 3190 each direction, with the results of the previous (and current) row being 3191 re-used. 3192 3193 That is after each row is 'Sync'ed' into the image, the next row makes use of 3194 those values as part of the calculation of the next row. It repeats, but 3195 going in the oppisite (bottom-up) direction. 3196 3197 Because of this 're-use of results' this function can not make use of multi- 3198 threaded, parellel processing. 3199*/ 3200static ssize_t MorphologyPrimitiveDirect(Image *image, 3201 const MorphologyMethod method,const KernelInfo *kernel, 3202 ExceptionInfo *exception) 3203{ 3204 CacheView 3205 *morphology_view, 3206 *image_view; 3207 3208 MagickBooleanType 3209 status; 3210 3211 MagickOffsetType 3212 progress; 3213 3214 OffsetInfo 3215 offset; 3216 3217 size_t 3218 width, 3219 changed; 3220 3221 ssize_t 3222 y; 3223 3224 assert(image != (Image *) NULL); 3225 assert(image->signature == MagickSignature); 3226 assert(kernel != (KernelInfo *) NULL); 3227 assert(kernel->signature == MagickSignature); 3228 assert(exception != (ExceptionInfo *) NULL); 3229 assert(exception->signature == MagickSignature); 3230 status=MagickTrue; 3231 changed=0; 3232 progress=0; 3233 switch(method) 3234 { 3235 case DistanceMorphology: 3236 case VoronoiMorphology: 3237 { 3238 /* 3239 Kernel reflected about origin. 3240 */ 3241 offset.x=(ssize_t) kernel->width-kernel->x-1; 3242 offset.y=(ssize_t) kernel->height-kernel->y-1; 3243 break; 3244 } 3245 default: 3246 { 3247 offset.x=kernel->x; 3248 offset.y=kernel->y; 3249 break; 3250 } 3251 } 3252 /* 3253 Two views into same image, do not thread. 3254 */ 3255 image_view=AcquireVirtualCacheView(image,exception); 3256 morphology_view=AcquireAuthenticCacheView(image,exception); 3257 width=image->columns+kernel->width-1; 3258 for (y=0; y < (ssize_t) image->rows; y++) 3259 { 3260 register const Quantum 3261 *restrict p; 3262 3263 register Quantum 3264 *restrict q; 3265 3266 register ssize_t 3267 x; 3268 3269 ssize_t 3270 center; 3271 3272 /* 3273 Read virtual pixels, and authentic pixels, from the same image! We read 3274 using virtual to get virtual pixel handling, but write back into the same 3275 image. 3276 3277 Only top half of kernel is processed as we do a single pass downward 3278 through the image iterating the distance function as we go. 3279 */ 3280 if (status == MagickFalse) 3281 continue; 3282 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t) 3283 offset.y+1,exception); 3284 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1, 3285 exception); 3286 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL)) 3287 { 3288 status=MagickFalse; 3289 continue; 3290 } 3291 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+ 3292 GetPixelChannels(image)*offset.x); 3293 for (x=0; x < (ssize_t) image->columns; x++) 3294 { 3295 register ssize_t 3296 i; 3297 3298 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 3299 { 3300 double 3301 pixel; 3302 3303 PixelTrait 3304 traits; 3305 3306 register const MagickRealType 3307 *restrict k; 3308 3309 register const Quantum 3310 *restrict pixels; 3311 3312 register ssize_t 3313 u; 3314 3315 ssize_t 3316 v; 3317 3318 traits=GetPixelChannelTraits(image,i); 3319 if (traits == UndefinedPixelTrait) 3320 continue; 3321 if (((traits & CopyPixelTrait) != 0) || 3322 (GetPixelReadMask(image,p+center) == 0)) 3323 continue; 3324 pixels=p; 3325 pixel=(double) QuantumRange; 3326 switch (method) 3327 { 3328 case DistanceMorphology: 3329 { 3330 k=(&kernel->values[kernel->width*kernel->height-1]); 3331 for (v=0; v <= offset.y; v++) 3332 { 3333 for (u=0; u < (ssize_t) kernel->width; u++) 3334 { 3335 if (IfNaN(*k) == MagickFalse) 3336 { 3337 if ((pixels[i]+(*k)) < pixel) 3338 pixel=(double) pixels[i]+(*k); 3339 } 3340 k--; 3341 pixels+=GetPixelChannels(image); 3342 } 3343 pixels+=(image->columns-1)*GetPixelChannels(image); 3344 } 3345 k=(&kernel->values[kernel->width*(kernel->y+1)-1]); 3346 pixels=q-offset.x*GetPixelChannels(image); 3347 for (u=0; u < offset.x; u++) 3348 { 3349 if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0)) 3350 { 3351 if ((pixels[i]+(*k)) < pixel) 3352 pixel=(double) pixels[i]+(*k); 3353 } 3354 k--; 3355 pixels+=GetPixelChannels(image); 3356 } 3357 break; 3358 } 3359 case VoronoiMorphology: 3360 { 3361 k=(&kernel->values[kernel->width*kernel->height-1]); 3362 for (v=0; v < offset.y; v++) 3363 { 3364 for (u=0; u < (ssize_t) kernel->width; u++) 3365 { 3366 if (IfNaN(*k) == MagickFalse) 3367 { 3368 if ((pixels[i]+(*k)) < pixel) 3369 pixel=(double) pixels[i]+(*k); 3370 } 3371 k--; 3372 pixels+=GetPixelChannels(image); 3373 } 3374 pixels+=(image->columns-1)*GetPixelChannels(image); 3375 } 3376 k=(&kernel->values[kernel->width*(kernel->y+1)-1]); 3377 pixels=q-offset.x*GetPixelChannels(image); 3378 for (u=0; u < offset.x; u++) 3379 { 3380 if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0)) 3381 { 3382 if ((pixels[i]+(*k)) < pixel) 3383 pixel=(double) pixels[i]+(*k); 3384 } 3385 k--; 3386 pixels+=GetPixelChannels(image); 3387 } 3388 break; 3389 } 3390 default: 3391 break; 3392 } 3393 if (fabs(pixel-q[i]) > MagickEpsilon) 3394 changed++; 3395 q[i]=ClampToQuantum(pixel); 3396 } 3397 p+=GetPixelChannels(image); 3398 q+=GetPixelChannels(image); 3399 } 3400 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse) 3401 status=MagickFalse; 3402 if (image->progress_monitor != (MagickProgressMonitor) NULL) 3403 { 3404 MagickBooleanType 3405 proceed; 3406 3407 proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows); 3408 if (proceed == MagickFalse) 3409 status=MagickFalse; 3410 } 3411 } 3412 morphology_view=DestroyCacheView(morphology_view); 3413 image_view=DestroyCacheView(image_view); 3414 /* 3415 Do the reverse pass through the image. 3416 */ 3417 image_view=AcquireVirtualCacheView(image,exception); 3418 morphology_view=AcquireAuthenticCacheView(image,exception); 3419 for (y=(ssize_t) image->rows-1; y >= 0; y--) 3420 { 3421 register const Quantum 3422 *restrict p; 3423 3424 register Quantum 3425 *restrict q; 3426 3427 register ssize_t 3428 x; 3429 3430 ssize_t 3431 center; 3432 3433 /* 3434 Read virtual pixels, and authentic pixels, from the same image. We 3435 read using virtual to get virtual pixel handling, but write back 3436 into the same image. 3437 3438 Only the bottom half of the kernel is processed as we up the image. 3439 */ 3440 if (status == MagickFalse) 3441 continue; 3442 p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t) 3443 kernel->y+1,exception); 3444 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1, 3445 exception); 3446 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL)) 3447 { 3448 status=MagickFalse; 3449 continue; 3450 } 3451 p+=(image->columns-1)*GetPixelChannels(image); 3452 q+=(image->columns-1)*GetPixelChannels(image); 3453 center=(ssize_t) (offset.x*GetPixelChannels(image)); 3454 for (x=(ssize_t) image->columns-1; x >= 0; x--) 3455 { 3456 register ssize_t 3457 i; 3458 3459 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 3460 { 3461 double 3462 pixel; 3463 3464 PixelTrait 3465 traits; 3466 3467 register const MagickRealType 3468 *restrict k; 3469 3470 register const Quantum 3471 *restrict pixels; 3472 3473 register ssize_t 3474 u; 3475 3476 ssize_t 3477 v; 3478 3479 traits=GetPixelChannelTraits(image,i); 3480 if (traits == UndefinedPixelTrait) 3481 continue; 3482 if (((traits & CopyPixelTrait) != 0) || 3483 (GetPixelReadMask(image,p+center) == 0)) 3484 continue; 3485 pixels=p; 3486 pixel=(double) QuantumRange; 3487 switch (method) 3488 { 3489 case DistanceMorphology: 3490 { 3491 k=(&kernel->values[kernel->width*(kernel->y+1)-1]); 3492 for (v=offset.y; v < (ssize_t) kernel->height; v++) 3493 { 3494 for (u=0; u < (ssize_t) kernel->width; u++) 3495 { 3496 if (IfNaN(*k) == MagickFalse) 3497 { 3498 if ((pixels[i]+(*k)) < pixel) 3499 pixel=(double) pixels[i]+(*k); 3500 } 3501 k--; 3502 pixels+=GetPixelChannels(image); 3503 } 3504 pixels+=(image->columns-1)*GetPixelChannels(image); 3505 } 3506 k=(&kernel->values[kernel->width*kernel->y+kernel->x-1]); 3507 pixels=q-offset.x*GetPixelChannels(image); 3508 for (u=offset.x+1; u < (ssize_t) kernel->width; u++) 3509 { 3510 if ((IfNaN(*k) == MagickFalse) && 3511 ((x+u-offset.x) < (ssize_t) image->columns)) 3512 { 3513 if ((pixels[i]+(*k)) < pixel) 3514 pixel=(double) pixels[i]+(*k); 3515 } 3516 k--; 3517 pixels+=GetPixelChannels(image); 3518 } 3519 break; 3520 } 3521 case VoronoiMorphology: 3522 { 3523 k=(&kernel->values[kernel->width*(kernel->y+1)-1]); 3524 for (v=offset.y; v < (ssize_t) kernel->height; v++) 3525 { 3526 for (u=0; u < (ssize_t) kernel->width; u++) 3527 { 3528 if (IfNaN(*k) == MagickFalse) 3529 { 3530 if ((pixels[i]+(*k)) < pixel) 3531 pixel=(double) pixels[i]+(*k); 3532 } 3533 k--; 3534 pixels+=GetPixelChannels(image); 3535 } 3536 pixels+=(image->columns-1)*GetPixelChannels(image); 3537 } 3538 k=(&kernel->values[kernel->width*(kernel->y+1)-1]); 3539 pixels=q-offset.x*GetPixelChannels(image); 3540 for (u=offset.x+1; u < (ssize_t) kernel->width; u++) 3541 { 3542 if ((IfNaN(*k) == MagickFalse) && 3543 ((x+u-offset.x) < (ssize_t) image->columns)) 3544 { 3545 if ((pixels[i]+(*k)) < pixel) 3546 pixel=(double) pixels[i]+(*k); 3547 } 3548 k--; 3549 pixels+=GetPixelChannels(image); 3550 } 3551 break; 3552 } 3553 default: 3554 break; 3555 } 3556 if (fabs(pixel-q[i]) > MagickEpsilon) 3557 changed++; 3558 q[i]=ClampToQuantum(pixel); 3559 } 3560 p-=GetPixelChannels(image); 3561 q-=GetPixelChannels(image); 3562 } 3563 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse) 3564 status=MagickFalse; 3565 if (image->progress_monitor != (MagickProgressMonitor) NULL) 3566 { 3567 MagickBooleanType 3568 proceed; 3569 3570 proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows); 3571 if (proceed == MagickFalse) 3572 status=MagickFalse; 3573 } 3574 } 3575 morphology_view=DestroyCacheView(morphology_view); 3576 image_view=DestroyCacheView(image_view); 3577 return(status ? (ssize_t) changed : -1); 3578} 3579 3580/* 3581 Apply a Morphology by calling one of the above low level primitive 3582 application functions. This function handles any iteration loops, 3583 composition or re-iteration of results, and compound morphology methods that 3584 is based on multiple low-level (staged) morphology methods. 3585 3586 Basically this provides the complex glue between the requested morphology 3587 method and raw low-level implementation (above). 3588*/ 3589MagickPrivate Image *MorphologyApply(const Image *image, 3590 const MorphologyMethod method, const ssize_t iterations, 3591 const KernelInfo *kernel, const CompositeOperator compose,const double bias, 3592 ExceptionInfo *exception) 3593{ 3594 CompositeOperator 3595 curr_compose; 3596 3597 Image 3598 *curr_image, /* Image we are working with or iterating */ 3599 *work_image, /* secondary image for primitive iteration */ 3600 *save_image, /* saved image - for 'edge' method only */ 3601 *rslt_image; /* resultant image - after multi-kernel handling */ 3602 3603 KernelInfo 3604 *reflected_kernel, /* A reflected copy of the kernel (if needed) */ 3605 *norm_kernel, /* the current normal un-reflected kernel */ 3606 *rflt_kernel, /* the current reflected kernel (if needed) */ 3607 *this_kernel; /* the kernel being applied */ 3608 3609 MorphologyMethod 3610 primitive; /* the current morphology primitive being applied */ 3611 3612 CompositeOperator 3613 rslt_compose; /* multi-kernel compose method for results to use */ 3614 3615 MagickBooleanType 3616 special, /* do we use a direct modify function? */ 3617 verbose; /* verbose output of results */ 3618 3619 size_t 3620 method_loop, /* Loop 1: number of compound method iterations (norm 1) */ 3621 method_limit, /* maximum number of compound method iterations */ 3622 kernel_number, /* Loop 2: the kernel number being applied */ 3623 stage_loop, /* Loop 3: primitive loop for compound morphology */ 3624 stage_limit, /* how many primitives are in this compound */ 3625 kernel_loop, /* Loop 4: iterate the kernel over image */ 3626 kernel_limit, /* number of times to iterate kernel */ 3627 count, /* total count of primitive steps applied */ 3628 kernel_changed, /* total count of changed using iterated kernel */ 3629 method_changed; /* total count of changed over method iteration */ 3630 3631 ssize_t 3632 changed; /* number pixels changed by last primitive operation */ 3633 3634 char 3635 v_info[80]; 3636 3637 assert(image != (Image *) NULL); 3638 assert(image->signature == MagickSignature); 3639 assert(kernel != (KernelInfo *) NULL); 3640 assert(kernel->signature == MagickSignature); 3641 assert(exception != (ExceptionInfo *) NULL); 3642 assert(exception->signature == MagickSignature); 3643 3644 count = 0; /* number of low-level morphology primitives performed */ 3645 if ( iterations == 0 ) 3646 return((Image *)NULL); /* null operation - nothing to do! */ 3647 3648 kernel_limit = (size_t) iterations; 3649 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */ 3650 kernel_limit = image->columns>image->rows ? image->columns : image->rows; 3651 3652 verbose = IsStringTrue(GetImageArtifact(image,"verbose")); 3653 3654 /* initialise for cleanup */ 3655 curr_image = (Image *) image; 3656 curr_compose = image->compose; 3657 (void) curr_compose; 3658 work_image = save_image = rslt_image = (Image *) NULL; 3659 reflected_kernel = (KernelInfo *) NULL; 3660 3661 /* Initialize specific methods 3662 * + which loop should use the given iteratations 3663 * + how many primitives make up the compound morphology 3664 * + multi-kernel compose method to use (by default) 3665 */ 3666 method_limit = 1; /* just do method once, unless otherwise set */ 3667 stage_limit = 1; /* assume method is not a compound */ 3668 special = MagickFalse; /* assume it is NOT a direct modify primitive */ 3669 rslt_compose = compose; /* and we are composing multi-kernels as given */ 3670 switch( method ) { 3671 case SmoothMorphology: /* 4 primitive compound morphology */ 3672 stage_limit = 4; 3673 break; 3674 case OpenMorphology: /* 2 primitive compound morphology */ 3675 case OpenIntensityMorphology: 3676 case TopHatMorphology: 3677 case CloseMorphology: 3678 case CloseIntensityMorphology: 3679 case BottomHatMorphology: 3680 case EdgeMorphology: 3681 stage_limit = 2; 3682 break; 3683 case HitAndMissMorphology: 3684 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */ 3685 /* FALL THUR */ 3686 case ThinningMorphology: 3687 case ThickenMorphology: 3688 method_limit = kernel_limit; /* iterate the whole method */ 3689 kernel_limit = 1; /* do not do kernel iteration */ 3690 break; 3691 case DistanceMorphology: 3692 case VoronoiMorphology: 3693 special = MagickTrue; /* use special direct primative */ 3694 break; 3695 default: 3696 break; 3697 } 3698 3699 /* Apply special methods with special requirments 3700 ** For example, single run only, or post-processing requirements 3701 */ 3702 if ( special == MagickTrue ) 3703 { 3704 rslt_image=CloneImage(image,0,0,MagickTrue,exception); 3705 if (rslt_image == (Image *) NULL) 3706 goto error_cleanup; 3707 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse) 3708 goto error_cleanup; 3709 3710 changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception); 3711 3712 if ( IfMagickTrue(verbose) ) 3713 (void) (void) FormatLocaleFile(stderr, 3714 "%s:%.20g.%.20g #%.20g => Changed %.20g\n", 3715 CommandOptionToMnemonic(MagickMorphologyOptions, method), 3716 1.0,0.0,1.0, (double) changed); 3717 3718 if ( changed < 0 ) 3719 goto error_cleanup; 3720 3721 if ( method == VoronoiMorphology ) { 3722 /* Preserve the alpha channel of input image - but turned it off */ 3723 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel, 3724 exception); 3725 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp, 3726 MagickTrue,0,0,exception); 3727 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel, 3728 exception); 3729 } 3730 goto exit_cleanup; 3731 } 3732 3733 /* Handle user (caller) specified multi-kernel composition method */ 3734 if ( compose != UndefinedCompositeOp ) 3735 rslt_compose = compose; /* override default composition for method */ 3736 if ( rslt_compose == UndefinedCompositeOp ) 3737 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */ 3738 3739 /* Some methods require a reflected kernel to use with primitives. 3740 * Create the reflected kernel for those methods. */ 3741 switch ( method ) { 3742 case CorrelateMorphology: 3743 case CloseMorphology: 3744 case CloseIntensityMorphology: 3745 case BottomHatMorphology: 3746 case SmoothMorphology: 3747 reflected_kernel = CloneKernelInfo(kernel); 3748 if (reflected_kernel == (KernelInfo *) NULL) 3749 goto error_cleanup; 3750 RotateKernelInfo(reflected_kernel,180); 3751 break; 3752 default: 3753 break; 3754 } 3755 3756 /* Loops around more primitive morpholgy methods 3757 ** erose, dilate, open, close, smooth, edge, etc... 3758 */ 3759 /* Loop 1: iterate the compound method */ 3760 method_loop = 0; 3761 method_changed = 1; 3762 while ( method_loop < method_limit && method_changed > 0 ) { 3763 method_loop++; 3764 method_changed = 0; 3765 3766 /* Loop 2: iterate over each kernel in a multi-kernel list */ 3767 norm_kernel = (KernelInfo *) kernel; 3768 this_kernel = (KernelInfo *) kernel; 3769 rflt_kernel = reflected_kernel; 3770 3771 kernel_number = 0; 3772 while ( norm_kernel != NULL ) { 3773 3774 /* Loop 3: Compound Morphology Staging - Select Primative to apply */ 3775 stage_loop = 0; /* the compound morphology stage number */ 3776 while ( stage_loop < stage_limit ) { 3777 stage_loop++; /* The stage of the compound morphology */ 3778 3779 /* Select primitive morphology for this stage of compound method */ 3780 this_kernel = norm_kernel; /* default use unreflected kernel */ 3781 primitive = method; /* Assume method is a primitive */ 3782 switch( method ) { 3783 case ErodeMorphology: /* just erode */ 3784 case EdgeInMorphology: /* erode and image difference */ 3785 primitive = ErodeMorphology; 3786 break; 3787 case DilateMorphology: /* just dilate */ 3788 case EdgeOutMorphology: /* dilate and image difference */ 3789 primitive = DilateMorphology; 3790 break; 3791 case OpenMorphology: /* erode then dialate */ 3792 case TopHatMorphology: /* open and image difference */ 3793 primitive = ErodeMorphology; 3794 if ( stage_loop == 2 ) 3795 primitive = DilateMorphology; 3796 break; 3797 case OpenIntensityMorphology: 3798 primitive = ErodeIntensityMorphology; 3799 if ( stage_loop == 2 ) 3800 primitive = DilateIntensityMorphology; 3801 break; 3802 case CloseMorphology: /* dilate, then erode */ 3803 case BottomHatMorphology: /* close and image difference */ 3804 this_kernel = rflt_kernel; /* use the reflected kernel */ 3805 primitive = DilateMorphology; 3806 if ( stage_loop == 2 ) 3807 primitive = ErodeMorphology; 3808 break; 3809 case CloseIntensityMorphology: 3810 this_kernel = rflt_kernel; /* use the reflected kernel */ 3811 primitive = DilateIntensityMorphology; 3812 if ( stage_loop == 2 ) 3813 primitive = ErodeIntensityMorphology; 3814 break; 3815 case SmoothMorphology: /* open, close */ 3816 switch ( stage_loop ) { 3817 case 1: /* start an open method, which starts with Erode */ 3818 primitive = ErodeMorphology; 3819 break; 3820 case 2: /* now Dilate the Erode */ 3821 primitive = DilateMorphology; 3822 break; 3823 case 3: /* Reflect kernel a close */ 3824 this_kernel = rflt_kernel; /* use the reflected kernel */ 3825 primitive = DilateMorphology; 3826 break; 3827 case 4: /* Finish the Close */ 3828 this_kernel = rflt_kernel; /* use the reflected kernel */ 3829 primitive = ErodeMorphology; 3830 break; 3831 } 3832 break; 3833 case EdgeMorphology: /* dilate and erode difference */ 3834 primitive = DilateMorphology; 3835 if ( stage_loop == 2 ) { 3836 save_image = curr_image; /* save the image difference */ 3837 curr_image = (Image *) image; 3838 primitive = ErodeMorphology; 3839 } 3840 break; 3841 case CorrelateMorphology: 3842 /* A Correlation is a Convolution with a reflected kernel. 3843 ** However a Convolution is a weighted sum using a reflected 3844 ** kernel. It may seem stange to convert a Correlation into a 3845 ** Convolution as the Correlation is the simplier method, but 3846 ** Convolution is much more commonly used, and it makes sense to 3847 ** implement it directly so as to avoid the need to duplicate the 3848 ** kernel when it is not required (which is typically the 3849 ** default). 3850 */ 3851 this_kernel = rflt_kernel; /* use the reflected kernel */ 3852 primitive = ConvolveMorphology; 3853 break; 3854 default: 3855 break; 3856 } 3857 assert( this_kernel != (KernelInfo *) NULL ); 3858 3859 /* Extra information for debugging compound operations */ 3860 if ( IfMagickTrue(verbose) ) { 3861 if ( stage_limit > 1 ) 3862 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ", 3863 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double) 3864 method_loop,(double) stage_loop); 3865 else if ( primitive != method ) 3866 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ", 3867 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double) 3868 method_loop); 3869 else 3870 v_info[0] = '\0'; 3871 } 3872 3873 /* Loop 4: Iterate the kernel with primitive */ 3874 kernel_loop = 0; 3875 kernel_changed = 0; 3876 changed = 1; 3877 while ( kernel_loop < kernel_limit && changed > 0 ) { 3878 kernel_loop++; /* the iteration of this kernel */ 3879 3880 /* Create a clone as the destination image, if not yet defined */ 3881 if ( work_image == (Image *) NULL ) 3882 { 3883 work_image=CloneImage(image,0,0,MagickTrue,exception); 3884 if (work_image == (Image *) NULL) 3885 goto error_cleanup; 3886 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse) 3887 goto error_cleanup; 3888 } 3889 3890 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */ 3891 count++; 3892 changed = MorphologyPrimitive(curr_image, work_image, primitive, 3893 this_kernel, bias, exception); 3894 3895 if ( IfMagickTrue(verbose) ) { 3896 if ( kernel_loop > 1 ) 3897 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */ 3898 (void) (void) FormatLocaleFile(stderr, 3899 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g", 3900 v_info,CommandOptionToMnemonic(MagickMorphologyOptions, 3901 primitive),(this_kernel == rflt_kernel ) ? "*" : "", 3902 (double) (method_loop+kernel_loop-1),(double) kernel_number, 3903 (double) count,(double) changed); 3904 } 3905 if ( changed < 0 ) 3906 goto error_cleanup; 3907 #pragma omp flush(changed) 3908 kernel_changed += changed; 3909 method_changed += changed; 3910 3911 /* prepare next loop */ 3912 { Image *tmp = work_image; /* swap images for iteration */ 3913 work_image = curr_image; 3914 curr_image = tmp; 3915 } 3916 if ( work_image == image ) 3917 work_image = (Image *) NULL; /* replace input 'image' */ 3918 3919 } /* End Loop 4: Iterate the kernel with primitive */ 3920 3921 if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed ) 3922 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed); 3923 if ( IfMagickTrue(verbose) && stage_loop < stage_limit ) 3924 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */ 3925 3926#if 0 3927 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image); 3928 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image); 3929 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image); 3930 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image); 3931 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image); 3932#endif 3933 3934 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */ 3935 3936 /* Final Post-processing for some Compound Methods 3937 ** 3938 ** The removal of any 'Sync' channel flag in the Image Compositon 3939 ** below ensures the methematical compose method is applied in a 3940 ** purely mathematical way, and only to the selected channels. 3941 ** Turn off SVG composition 'alpha blending'. 3942 */ 3943 switch( method ) { 3944 case EdgeOutMorphology: 3945 case EdgeInMorphology: 3946 case TopHatMorphology: 3947 case BottomHatMorphology: 3948 if ( IfMagickTrue(verbose) ) 3949 (void) FormatLocaleFile(stderr, 3950 "\n%s: Difference with original image",CommandOptionToMnemonic( 3951 MagickMorphologyOptions, method) ); 3952 (void) CompositeImage(curr_image,image,DifferenceCompositeOp, 3953 MagickTrue,0,0,exception); 3954 break; 3955 case EdgeMorphology: 3956 if ( IfMagickTrue(verbose) ) 3957 (void) FormatLocaleFile(stderr, 3958 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic( 3959 MagickMorphologyOptions, method) ); 3960 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp, 3961 MagickTrue,0,0,exception); 3962 save_image = DestroyImage(save_image); /* finished with save image */ 3963 break; 3964 default: 3965 break; 3966 } 3967 3968 /* multi-kernel handling: re-iterate, or compose results */ 3969 if ( kernel->next == (KernelInfo *) NULL ) 3970 rslt_image = curr_image; /* just return the resulting image */ 3971 else if ( rslt_compose == NoCompositeOp ) 3972 { if ( IfMagickTrue(verbose) ) { 3973 if ( this_kernel->next != (KernelInfo *) NULL ) 3974 (void) FormatLocaleFile(stderr, " (re-iterate)"); 3975 else 3976 (void) FormatLocaleFile(stderr, " (done)"); 3977 } 3978 rslt_image = curr_image; /* return result, and re-iterate */ 3979 } 3980 else if ( rslt_image == (Image *) NULL) 3981 { if ( IfMagickTrue(verbose) ) 3982 (void) FormatLocaleFile(stderr, " (save for compose)"); 3983 rslt_image = curr_image; 3984 curr_image = (Image *) image; /* continue with original image */ 3985 } 3986 else 3987 { /* Add the new 'current' result to the composition 3988 ** 3989 ** The removal of any 'Sync' channel flag in the Image Compositon 3990 ** below ensures the methematical compose method is applied in a 3991 ** purely mathematical way, and only to the selected channels. 3992 ** IE: Turn off SVG composition 'alpha blending'. 3993 */ 3994 if ( IfMagickTrue(verbose) ) 3995 (void) FormatLocaleFile(stderr, " (compose \"%s\")", 3996 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) ); 3997 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue, 3998 0,0,exception); 3999 curr_image = DestroyImage(curr_image); 4000 curr_image = (Image *) image; /* continue with original image */ 4001 } 4002 if ( IfMagickTrue(verbose) ) 4003 (void) FormatLocaleFile(stderr, "\n"); 4004 4005 /* loop to the next kernel in a multi-kernel list */ 4006 norm_kernel = norm_kernel->next; 4007 if ( rflt_kernel != (KernelInfo *) NULL ) 4008 rflt_kernel = rflt_kernel->next; 4009 kernel_number++; 4010 } /* End Loop 2: Loop over each kernel */ 4011 4012 } /* End Loop 1: compound method interation */ 4013 4014 goto exit_cleanup; 4015 4016 /* Yes goto's are bad, but it makes cleanup lot more efficient */ 4017error_cleanup: 4018 if ( curr_image == rslt_image ) 4019 curr_image = (Image *) NULL; 4020 if ( rslt_image != (Image *) NULL ) 4021 rslt_image = DestroyImage(rslt_image); 4022exit_cleanup: 4023 if ( curr_image == rslt_image || curr_image == image ) 4024 curr_image = (Image *) NULL; 4025 if ( curr_image != (Image *) NULL ) 4026 curr_image = DestroyImage(curr_image); 4027 if ( work_image != (Image *) NULL ) 4028 work_image = DestroyImage(work_image); 4029 if ( save_image != (Image *) NULL ) 4030 save_image = DestroyImage(save_image); 4031 if ( reflected_kernel != (KernelInfo *) NULL ) 4032 reflected_kernel = DestroyKernelInfo(reflected_kernel); 4033 return(rslt_image); 4034} 4035 4036 4037/* 4038%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4039% % 4040% % 4041% % 4042% M o r p h o l o g y I m a g e % 4043% % 4044% % 4045% % 4046%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4047% 4048% MorphologyImage() applies a user supplied kernel to the image according to 4049% the given mophology method. 4050% 4051% This function applies any and all user defined settings before calling 4052% the above internal function MorphologyApply(). 4053% 4054% User defined settings include... 4055% * Output Bias for Convolution and correlation ("-define convolve:bias=??") 4056% * Kernel Scale/normalize settings ("-define convolve:scale=??") 4057% This can also includes the addition of a scaled unity kernel. 4058% * Show Kernel being applied ("-define showkernel=1") 4059% 4060% Other operators that do not want user supplied options interfering, 4061% especially "convolve:bias" and "showkernel" should use MorphologyApply() 4062% directly. 4063% 4064% The format of the MorphologyImage method is: 4065% 4066% Image *MorphologyImage(const Image *image,MorphologyMethod method, 4067% const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception) 4068% 4069% A description of each parameter follows: 4070% 4071% o image: the image. 4072% 4073% o method: the morphology method to be applied. 4074% 4075% o iterations: apply the operation this many times (or no change). 4076% A value of -1 means loop until no change found. 4077% How this is applied may depend on the morphology method. 4078% Typically this is a value of 1. 4079% 4080% o kernel: An array of double representing the morphology kernel. 4081% Warning: kernel may be normalized for the Convolve method. 4082% 4083% o exception: return any errors or warnings in this structure. 4084% 4085*/ 4086MagickExport Image *MorphologyImage(const Image *image, 4087 const MorphologyMethod method,const ssize_t iterations, 4088 const KernelInfo *kernel,ExceptionInfo *exception) 4089{ 4090 KernelInfo 4091 *curr_kernel; 4092 4093 CompositeOperator 4094 compose; 4095 4096 Image 4097 *morphology_image; 4098 4099 double 4100 bias; 4101 4102 curr_kernel = (KernelInfo *) kernel; 4103 bias=0.0; 4104 compose = UndefinedCompositeOp; /* use default for method */ 4105 4106 /* Apply Convolve/Correlate Normalization and Scaling Factors. 4107 * This is done BEFORE the ShowKernelInfo() function is called so that 4108 * users can see the results of the 'option:convolve:scale' option. 4109 */ 4110 if ( method == ConvolveMorphology || method == CorrelateMorphology ) { 4111 const char 4112 *artifact; 4113 4114 /* Get the bias value as it will be needed */ 4115 artifact = GetImageArtifact(image,"convolve:bias"); 4116 if ( artifact != (const char *) NULL) { 4117 if (IfMagickFalse(IsGeometry(artifact))) 4118 (void) ThrowMagickException(exception,GetMagickModule(), 4119 OptionWarning,"InvalidSetting","'%s' '%s'", 4120 "convolve:bias",artifact); 4121 else 4122 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0); 4123 } 4124 4125 /* Scale kernel according to user wishes */ 4126 artifact = GetImageArtifact(image,"convolve:scale"); 4127 if ( artifact != (const char *)NULL ) { 4128 if (IfMagickFalse(IsGeometry(artifact))) 4129 (void) ThrowMagickException(exception,GetMagickModule(), 4130 OptionWarning,"InvalidSetting","'%s' '%s'", 4131 "convolve:scale",artifact); 4132 else { 4133 if ( curr_kernel == kernel ) 4134 curr_kernel = CloneKernelInfo(kernel); 4135 if (curr_kernel == (KernelInfo *) NULL) 4136 return((Image *) NULL); 4137 ScaleGeometryKernelInfo(curr_kernel, artifact); 4138 } 4139 } 4140 } 4141 4142 /* display the (normalized) kernel via stderr */ 4143 if ( IfStringTrue(GetImageArtifact(image,"showkernel")) 4144 || IfStringTrue(GetImageArtifact(image,"convolve:showkernel")) 4145 || IfStringTrue(GetImageArtifact(image,"morphology:showkernel")) ) 4146 ShowKernelInfo(curr_kernel); 4147 4148 /* Override the default handling of multi-kernel morphology results 4149 * If 'Undefined' use the default method 4150 * If 'None' (default for 'Convolve') re-iterate previous result 4151 * Otherwise merge resulting images using compose method given. 4152 * Default for 'HitAndMiss' is 'Lighten'. 4153 */ 4154 { const char 4155 *artifact; 4156 ssize_t 4157 parse; 4158 4159 artifact = GetImageArtifact(image,"morphology:compose"); 4160 if ( artifact != (const char *) NULL) { 4161 parse=ParseCommandOption(MagickComposeOptions, 4162 MagickFalse,artifact); 4163 if ( parse < 0 ) 4164 (void) ThrowMagickException(exception,GetMagickModule(), 4165 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'", 4166 "morphology:compose",artifact); 4167 else 4168 compose=(CompositeOperator)parse; 4169 } 4170 } 4171 /* Apply the Morphology */ 4172 morphology_image = MorphologyApply(image,method,iterations, 4173 curr_kernel,compose,bias,exception); 4174 4175 /* Cleanup and Exit */ 4176 if ( curr_kernel != kernel ) 4177 curr_kernel=DestroyKernelInfo(curr_kernel); 4178 return(morphology_image); 4179} 4180 4181/* 4182%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4183% % 4184% % 4185% % 4186+ R o t a t e K e r n e l I n f o % 4187% % 4188% % 4189% % 4190%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4191% 4192% RotateKernelInfo() rotates the kernel by the angle given. 4193% 4194% Currently it is restricted to 90 degree angles, of either 1D kernels 4195% or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels. 4196% It will ignore usless rotations for specific 'named' built-in kernels. 4197% 4198% The format of the RotateKernelInfo method is: 4199% 4200% void RotateKernelInfo(KernelInfo *kernel, double angle) 4201% 4202% A description of each parameter follows: 4203% 4204% o kernel: the Morphology/Convolution kernel 4205% 4206% o angle: angle to rotate in degrees 4207% 4208% This function is currently internal to this module only, but can be exported 4209% to other modules if needed. 4210*/ 4211static void RotateKernelInfo(KernelInfo *kernel, double angle) 4212{ 4213 /* angle the lower kernels first */ 4214 if ( kernel->next != (KernelInfo *) NULL) 4215 RotateKernelInfo(kernel->next, angle); 4216 4217 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical 4218 ** 4219 ** TODO: expand beyond simple 90 degree rotates, flips and flops 4220 */ 4221 4222 /* Modulus the angle */ 4223 angle = fmod(angle, 360.0); 4224 if ( angle < 0 ) 4225 angle += 360.0; 4226 4227 if ( 337.5 < angle || angle <= 22.5 ) 4228 return; /* Near zero angle - no change! - At least not at this time */ 4229 4230 /* Handle special cases */ 4231 switch (kernel->type) { 4232 /* These built-in kernels are cylindrical kernels, rotating is useless */ 4233 case GaussianKernel: 4234 case DoGKernel: 4235 case LoGKernel: 4236 case DiskKernel: 4237 case PeaksKernel: 4238 case LaplacianKernel: 4239 case ChebyshevKernel: 4240 case ManhattanKernel: 4241 case EuclideanKernel: 4242 return; 4243 4244 /* These may be rotatable at non-90 angles in the future */ 4245 /* but simply rotating them in multiples of 90 degrees is useless */ 4246 case SquareKernel: 4247 case DiamondKernel: 4248 case PlusKernel: 4249 case CrossKernel: 4250 return; 4251 4252 /* These only allows a +/-90 degree rotation (by transpose) */ 4253 /* A 180 degree rotation is useless */ 4254 case BlurKernel: 4255 if ( 135.0 < angle && angle <= 225.0 ) 4256 return; 4257 if ( 225.0 < angle && angle <= 315.0 ) 4258 angle -= 180; 4259 break; 4260 4261 default: 4262 break; 4263 } 4264 /* Attempt rotations by 45 degrees -- 3x3 kernels only */ 4265 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 ) 4266 { 4267 if ( kernel->width == 3 && kernel->height == 3 ) 4268 { /* Rotate a 3x3 square by 45 degree angle */ 4269 double t = kernel->values[0]; 4270 kernel->values[0] = kernel->values[3]; 4271 kernel->values[3] = kernel->values[6]; 4272 kernel->values[6] = kernel->values[7]; 4273 kernel->values[7] = kernel->values[8]; 4274 kernel->values[8] = kernel->values[5]; 4275 kernel->values[5] = kernel->values[2]; 4276 kernel->values[2] = kernel->values[1]; 4277 kernel->values[1] = t; 4278 /* rotate non-centered origin */ 4279 if ( kernel->x != 1 || kernel->y != 1 ) { 4280 ssize_t x,y; 4281 x = (ssize_t) kernel->x-1; 4282 y = (ssize_t) kernel->y-1; 4283 if ( x == y ) x = 0; 4284 else if ( x == 0 ) x = -y; 4285 else if ( x == -y ) y = 0; 4286 else if ( y == 0 ) y = x; 4287 kernel->x = (ssize_t) x+1; 4288 kernel->y = (ssize_t) y+1; 4289 } 4290 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */ 4291 kernel->angle = fmod(kernel->angle+45.0, 360.0); 4292 } 4293 else 4294 perror("Unable to rotate non-3x3 kernel by 45 degrees"); 4295 } 4296 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 ) 4297 { 4298 if ( kernel->width == 1 || kernel->height == 1 ) 4299 { /* Do a transpose of a 1 dimensional kernel, 4300 ** which results in a fast 90 degree rotation of some type. 4301 */ 4302 ssize_t 4303 t; 4304 t = (ssize_t) kernel->width; 4305 kernel->width = kernel->height; 4306 kernel->height = (size_t) t; 4307 t = kernel->x; 4308 kernel->x = kernel->y; 4309 kernel->y = t; 4310 if ( kernel->width == 1 ) { 4311 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */ 4312 kernel->angle = fmod(kernel->angle+90.0, 360.0); 4313 } else { 4314 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */ 4315 kernel->angle = fmod(kernel->angle+270.0, 360.0); 4316 } 4317 } 4318 else if ( kernel->width == kernel->height ) 4319 { /* Rotate a square array of values by 90 degrees */ 4320 { register ssize_t 4321 i,j,x,y; 4322 4323 register MagickRealType 4324 *k,t; 4325 4326 k=kernel->values; 4327 for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--) 4328 for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--) 4329 { t = k[i+j*kernel->width]; 4330 k[i+j*kernel->width] = k[j+x*kernel->width]; 4331 k[j+x*kernel->width] = k[x+y*kernel->width]; 4332 k[x+y*kernel->width] = k[y+i*kernel->width]; 4333 k[y+i*kernel->width] = t; 4334 } 4335 } 4336 /* rotate the origin - relative to center of array */ 4337 { register ssize_t x,y; 4338 x = (ssize_t) (kernel->x*2-kernel->width+1); 4339 y = (ssize_t) (kernel->y*2-kernel->height+1); 4340 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2; 4341 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2; 4342 } 4343 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */ 4344 kernel->angle = fmod(kernel->angle+90.0, 360.0); 4345 } 4346 else 4347 perror("Unable to rotate a non-square, non-linear kernel 90 degrees"); 4348 } 4349 if ( 135.0 < angle && angle <= 225.0 ) 4350 { 4351 /* For a 180 degree rotation - also know as a reflection 4352 * This is actually a very very common operation! 4353 * Basically all that is needed is a reversal of the kernel data! 4354 * And a reflection of the origon 4355 */ 4356 MagickRealType 4357 t; 4358 4359 register MagickRealType 4360 *k; 4361 4362 ssize_t 4363 i, 4364 j; 4365 4366 k=kernel->values; 4367 j=(ssize_t) (kernel->width*kernel->height-1); 4368 for (i=0; i < j; i++, j--) 4369 t=k[i], k[i]=k[j], k[j]=t; 4370 4371 kernel->x = (ssize_t) kernel->width - kernel->x - 1; 4372 kernel->y = (ssize_t) kernel->height - kernel->y - 1; 4373 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */ 4374 kernel->angle = fmod(kernel->angle+180.0, 360.0); 4375 } 4376 /* At this point angle should at least between -45 (315) and +45 degrees 4377 * In the future some form of non-orthogonal angled rotates could be 4378 * performed here, posibily with a linear kernel restriction. 4379 */ 4380 4381 return; 4382} 4383 4384/* 4385%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4386% % 4387% % 4388% % 4389% S c a l e G e o m e t r y K e r n e l I n f o % 4390% % 4391% % 4392% % 4393%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4394% 4395% ScaleGeometryKernelInfo() takes a geometry argument string, typically 4396% provided as a "-set option:convolve:scale {geometry}" user setting, 4397% and modifies the kernel according to the parsed arguments of that setting. 4398% 4399% The first argument (and any normalization flags) are passed to 4400% ScaleKernelInfo() to scale/normalize the kernel. The second argument 4401% is then passed to UnityAddKernelInfo() to add a scled unity kernel 4402% into the scaled/normalized kernel. 4403% 4404% The format of the ScaleGeometryKernelInfo method is: 4405% 4406% void ScaleGeometryKernelInfo(KernelInfo *kernel, 4407% const double scaling_factor,const MagickStatusType normalize_flags) 4408% 4409% A description of each parameter follows: 4410% 4411% o kernel: the Morphology/Convolution kernel to modify 4412% 4413% o geometry: 4414% The geometry string to parse, typically from the user provided 4415% "-set option:convolve:scale {geometry}" setting. 4416% 4417*/ 4418MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel, 4419 const char *geometry) 4420{ 4421 MagickStatusType 4422 flags; 4423 4424 GeometryInfo 4425 args; 4426 4427 SetGeometryInfo(&args); 4428 flags = ParseGeometry(geometry, &args); 4429 4430#if 0 4431 /* For Debugging Geometry Input */ 4432 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n", 4433 flags, args.rho, args.sigma, args.xi, args.psi ); 4434#endif 4435 4436 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/ 4437 args.rho *= 0.01, args.sigma *= 0.01; 4438 4439 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */ 4440 args.rho = 1.0; 4441 if ( (flags & SigmaValue) == 0 ) 4442 args.sigma = 0.0; 4443 4444 /* Scale/Normalize the input kernel */ 4445 ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags); 4446 4447 /* Add Unity Kernel, for blending with original */ 4448 if ( (flags & SigmaValue) != 0 ) 4449 UnityAddKernelInfo(kernel, args.sigma); 4450 4451 return; 4452} 4453/* 4454%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4455% % 4456% % 4457% % 4458% S c a l e K e r n e l I n f o % 4459% % 4460% % 4461% % 4462%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4463% 4464% ScaleKernelInfo() scales the given kernel list by the given amount, with or 4465% without normalization of the sum of the kernel values (as per given flags). 4466% 4467% By default (no flags given) the values within the kernel is scaled 4468% directly using given scaling factor without change. 4469% 4470% If either of the two 'normalize_flags' are given the kernel will first be 4471% normalized and then further scaled by the scaling factor value given. 4472% 4473% Kernel normalization ('normalize_flags' given) is designed to ensure that 4474% any use of the kernel scaling factor with 'Convolve' or 'Correlate' 4475% morphology methods will fall into -1.0 to +1.0 range. Note that for 4476% non-HDRI versions of IM this may cause images to have any negative results 4477% clipped, unless some 'bias' is used. 4478% 4479% More specifically. Kernels which only contain positive values (such as a 4480% 'Gaussian' kernel) will be scaled so that those values sum to +1.0, 4481% ensuring a 0.0 to +1.0 output range for non-HDRI images. 4482% 4483% For Kernels that contain some negative values, (such as 'Sharpen' kernels) 4484% the kernel will be scaled by the absolute of the sum of kernel values, so 4485% that it will generally fall within the +/- 1.0 range. 4486% 4487% For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel 4488% will be scaled by just the sum of the postive values, so that its output 4489% range will again fall into the +/- 1.0 range. 4490% 4491% For special kernels designed for locating shapes using 'Correlate', (often 4492% only containing +1 and -1 values, representing foreground/brackground 4493% matching) a special normalization method is provided to scale the positive 4494% values separately to those of the negative values, so the kernel will be 4495% forced to become a zero-sum kernel better suited to such searches. 4496% 4497% WARNING: Correct normalization of the kernel assumes that the '*_range' 4498% attributes within the kernel structure have been correctly set during the 4499% kernels creation. 4500% 4501% NOTE: The values used for 'normalize_flags' have been selected specifically 4502% to match the use of geometry options, so that '!' means NormalizeValue, '^' 4503% means CorrelateNormalizeValue. All other GeometryFlags values are ignored. 4504% 4505% The format of the ScaleKernelInfo method is: 4506% 4507% void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor, 4508% const MagickStatusType normalize_flags ) 4509% 4510% A description of each parameter follows: 4511% 4512% o kernel: the Morphology/Convolution kernel 4513% 4514% o scaling_factor: 4515% multiply all values (after normalization) by this factor if not 4516% zero. If the kernel is normalized regardless of any flags. 4517% 4518% o normalize_flags: 4519% GeometryFlags defining normalization method to use. 4520% specifically: NormalizeValue, CorrelateNormalizeValue, 4521% and/or PercentValue 4522% 4523*/ 4524MagickExport void ScaleKernelInfo(KernelInfo *kernel, 4525 const double scaling_factor,const GeometryFlags normalize_flags) 4526{ 4527 register double 4528 pos_scale, 4529 neg_scale; 4530 4531 register ssize_t 4532 i; 4533 4534 /* do the other kernels in a multi-kernel list first */ 4535 if ( kernel->next != (KernelInfo *) NULL) 4536 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags); 4537 4538 /* Normalization of Kernel */ 4539 pos_scale = 1.0; 4540 if ( (normalize_flags&NormalizeValue) != 0 ) { 4541 if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon ) 4542 /* non-zero-summing kernel (generally positive) */ 4543 pos_scale = fabs(kernel->positive_range + kernel->negative_range); 4544 else 4545 /* zero-summing kernel */ 4546 pos_scale = kernel->positive_range; 4547 } 4548 /* Force kernel into a normalized zero-summing kernel */ 4549 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) { 4550 pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon ) 4551 ? kernel->positive_range : 1.0; 4552 neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon ) 4553 ? -kernel->negative_range : 1.0; 4554 } 4555 else 4556 neg_scale = pos_scale; 4557 4558 /* finialize scaling_factor for positive and negative components */ 4559 pos_scale = scaling_factor/pos_scale; 4560 neg_scale = scaling_factor/neg_scale; 4561 4562 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++) 4563 if ( ! IfNaN(kernel->values[i]) ) 4564 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale; 4565 4566 /* convolution output range */ 4567 kernel->positive_range *= pos_scale; 4568 kernel->negative_range *= neg_scale; 4569 /* maximum and minimum values in kernel */ 4570 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale; 4571 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale; 4572 4573 /* swap kernel settings if user's scaling factor is negative */ 4574 if ( scaling_factor < MagickEpsilon ) { 4575 double t; 4576 t = kernel->positive_range; 4577 kernel->positive_range = kernel->negative_range; 4578 kernel->negative_range = t; 4579 t = kernel->maximum; 4580 kernel->maximum = kernel->minimum; 4581 kernel->minimum = 1; 4582 } 4583 4584 return; 4585} 4586 4587/* 4588%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4589% % 4590% % 4591% % 4592% S h o w K e r n e l I n f o % 4593% % 4594% % 4595% % 4596%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4597% 4598% ShowKernelInfo() outputs the details of the given kernel defination to 4599% standard error, generally due to a users 'showkernel' option request. 4600% 4601% The format of the ShowKernel method is: 4602% 4603% void ShowKernelInfo(const KernelInfo *kernel) 4604% 4605% A description of each parameter follows: 4606% 4607% o kernel: the Morphology/Convolution kernel 4608% 4609*/ 4610MagickPrivate void ShowKernelInfo(const KernelInfo *kernel) 4611{ 4612 const KernelInfo 4613 *k; 4614 4615 size_t 4616 c, i, u, v; 4617 4618 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) { 4619 4620 (void) FormatLocaleFile(stderr, "Kernel"); 4621 if ( kernel->next != (KernelInfo *) NULL ) 4622 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c ); 4623 (void) FormatLocaleFile(stderr, " \"%s", 4624 CommandOptionToMnemonic(MagickKernelOptions, k->type) ); 4625 if ( fabs(k->angle) >= MagickEpsilon ) 4626 (void) FormatLocaleFile(stderr, "@%lg", k->angle); 4627 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long) 4628 k->width,(unsigned long) k->height,(long) k->x,(long) k->y); 4629 (void) FormatLocaleFile(stderr, 4630 " with values from %.*lg to %.*lg\n", 4631 GetMagickPrecision(), k->minimum, 4632 GetMagickPrecision(), k->maximum); 4633 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg", 4634 GetMagickPrecision(), k->negative_range, 4635 GetMagickPrecision(), k->positive_range); 4636 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon ) 4637 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n"); 4638 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon ) 4639 (void) FormatLocaleFile(stderr, " (Normalized)\n"); 4640 else 4641 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n", 4642 GetMagickPrecision(), k->positive_range+k->negative_range); 4643 for (i=v=0; v < k->height; v++) { 4644 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v ); 4645 for (u=0; u < k->width; u++, i++) 4646 if ( IfNaN(k->values[i]) ) 4647 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan"); 4648 else 4649 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3, 4650 GetMagickPrecision(), (double) k->values[i]); 4651 (void) FormatLocaleFile(stderr,"\n"); 4652 } 4653 } 4654} 4655 4656/* 4657%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4658% % 4659% % 4660% % 4661% U n i t y A d d K e r n a l I n f o % 4662% % 4663% % 4664% % 4665%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4666% 4667% UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel 4668% to the given pre-scaled and normalized Kernel. This in effect adds that 4669% amount of the original image into the resulting convolution kernel. This 4670% value is usually provided by the user as a percentage value in the 4671% 'convolve:scale' setting. 4672% 4673% The resulting effect is to convert the defined kernels into blended 4674% soft-blurs, unsharp kernels or into sharpening kernels. 4675% 4676% The format of the UnityAdditionKernelInfo method is: 4677% 4678% void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale ) 4679% 4680% A description of each parameter follows: 4681% 4682% o kernel: the Morphology/Convolution kernel 4683% 4684% o scale: 4685% scaling factor for the unity kernel to be added to 4686% the given kernel. 4687% 4688*/ 4689MagickExport void UnityAddKernelInfo(KernelInfo *kernel, 4690 const double scale) 4691{ 4692 /* do the other kernels in a multi-kernel list first */ 4693 if ( kernel->next != (KernelInfo *) NULL) 4694 UnityAddKernelInfo(kernel->next, scale); 4695 4696 /* Add the scaled unity kernel to the existing kernel */ 4697 kernel->values[kernel->x+kernel->y*kernel->width] += scale; 4698 CalcKernelMetaData(kernel); /* recalculate the meta-data */ 4699 4700 return; 4701} 4702 4703/* 4704%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4705% % 4706% % 4707% % 4708% Z e r o K e r n e l N a n s % 4709% % 4710% % 4711% % 4712%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4713% 4714% ZeroKernelNans() replaces any special 'nan' value that may be present in 4715% the kernel with a zero value. This is typically done when the kernel will 4716% be used in special hardware (GPU) convolution processors, to simply 4717% matters. 4718% 4719% The format of the ZeroKernelNans method is: 4720% 4721% void ZeroKernelNans (KernelInfo *kernel) 4722% 4723% A description of each parameter follows: 4724% 4725% o kernel: the Morphology/Convolution kernel 4726% 4727*/ 4728MagickPrivate void ZeroKernelNans(KernelInfo *kernel) 4729{ 4730 register size_t 4731 i; 4732 4733 /* do the other kernels in a multi-kernel list first */ 4734 if ( kernel->next != (KernelInfo *) NULL) 4735 ZeroKernelNans(kernel->next); 4736 4737 for (i=0; i < (kernel->width*kernel->height); i++) 4738 if ( IfNaN(kernel->values[i]) ) 4739 kernel->values[i] = 0.0; 4740 4741 return; 4742} 4743