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