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