« back to all changes in this revision
Viewing changes to magick/morphology.c
- browse files at revision 29
- compare with another revision
- download diff
- download tarball
- view history from revision 29
- Committer: Bazaar Package Importer
- Author(s): Colin Watson
- Date: 2011-06-15 11:05:28 UTC
- mfrom: (6.2.11 sid)
- Revision ID: james.westby@ubuntu.com-20110615110528-08jgo07a4846xh8d
Tags: 8:6.6.0.4-3ubuntu1
* Resynchronise with Debian (LP: #797595). Remaining changes:
- Make ufraw-batch (universe) a suggestion instead of a recommendation.
- Make debian/rules install target depend on check; they cannot reliably
be run in parallel.
- Don't set MAKEFLAGS in debian/rules; just pass it to the build.
- Make ufraw-batch (universe) a suggestion instead of a recommendation.
- Make debian/rules install target depend on check; they cannot reliably
be run in parallel.
- Don't set MAKEFLAGS in debian/rules; just pass it to the build.
- files added:
- m4/ac_compile_warnings.m4
- files removed:
- coders/debug.c
- files modified:
- .pc/601824.patch/magick/configure.c
- magick.sh.in *
added
removed
magick/morphology.c
33
33
% %
34
34
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35
35
%
36
% Morpology is the the application of various kernels, of any size and even
36
% Morpology is the the application of various kernals, of any size and even
37
37
% shape, to a image in various ways (typically binary, but not always).
38
38
%
39
39
% Convolution (weighted sum or average) is just one specific type of
65
65
#include "magick/memory_.h"
66
66
#include "magick/monitor-private.h"
67
67
#include "magick/morphology.h"
68
#include "magick/morphology-private.h"
69
68
#include "magick/option.h"
70
69
#include "magick/pixel-private.h"
71
70
#include "magick/prepress.h"
78
77
#include "magick/string-private.h"
79
78
#include "magick/token.h"
80
79
81
82
80
/*
83
** The following test is for special floating point numbers of value NaN (not
84
** a number), that may be used within a Kernel Definition. NaN's are defined
85
** as part of the IEEE standard for floating point number representation.
86
**
87
** These are used as a Kernel value to mean that this kernel position is not
88
** part of the kernel neighbourhood for convolution or morphology processing,
89
** and thus should be ignored. This allows the use of 'shaped' kernels.
90
**
91
** The special properity that two NaN's are never equal, even if they are from
92
** the same variable allow you to test if a value is special NaN value.
93
**
94
** This macro IsNaN() is thus is only true if the value given is NaN.
81
The following test is for special floating point numbers of value NaN (not
82
a number), that may be used within a Kernel Definition. NaN's are defined
83
as part of the IEEE standard for floating point number representation.
84
85
These are used a Kernel value of NaN means that that kernal position is not
86
part of the normal convolution or morphology process, and thus allowing the
87
use of 'shaped' kernels.
88
89
Special properities two NaN's are never equal, even if they are from the
90
same variable That is the IsNaN() macro is only true if the value is NaN.
95
91
*/
96
92
#define IsNan(a) ((a)!=(a))
97
93
111
107
112
108
/* Currently these are only internal to this module */
113
109
static void
114
CalcKernelMetaData(KernelInfo *),
115
ExpandMirrorKernelInfo(KernelInfo *),
116
ExpandRotateKernelInfo(KernelInfo *, const double),
117
110
RotateKernelInfo(KernelInfo *, double);
118
111
119
120
/* Quick function to find last kernel in a kernel list */
121
static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
122
{
123
while (kernel->next != (KernelInfo *) NULL)
124
kernel = kernel->next;
125
return(kernel);
126
}
127
128
129
112
/*
130
113
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131
114
% %
148
131
% anywhere within that array of values.
149
132
%
150
133
% Previously IM was restricted to a square of odd size using the exact
151
% center as origin, this is no ssize_ter the case, and any rectangular kernel
134
% center as origin, this is no longer the case, and any rectangular kernel
152
135
% with any value being declared the origin. This in turn allows the use of
153
136
% highly asymmetrical kernels.
154
137
%
159
142
% shape. However at least one non-nan value must be provided for correct
160
143
% working of a kernel.
161
144
%
162
% The returned kernel should be freed using the DestroyKernelInfo() when you
163
% are finished with it. Do not free this memory yourself.
145
% The returned kernel should be free using the DestroyKernelInfo() when you
146
% are finished with it.
164
147
%
165
148
% Input kernel defintion strings can consist of any of three types.
166
149
%
167
% "name:args[[@><]"
150
% "name:args"
168
151
% Select from one of the built in kernels, using the name and
169
152
% geometry arguments supplied. See AcquireKernelBuiltIn()
170
153
%
171
% "WxH[+X+Y][@><]:num, num, num ..."
172
% a kernel of size W by H, with W*H floating point numbers following.
154
% "WxH[+X+Y]:num, num, num ..."
155
% a kernal of size W by H, with W*H floating point numbers following.
173
156
% the 'center' can be optionally be defined at +X+Y (such that +0+0
174
157
% is top left corner). If not defined the pixel in the center, for
175
158
% odd sizes, or to the immediate top or left of center for even sizes
181
164
% square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
182
165
% Values can be space or comma separated. This is not recommended.
183
166
%
184
% You can define a 'list of kernels' which can be used by some morphology
185
% operators A list is defined as a semi-colon seperated list kernels.
186
%
187
% " kernel ; kernel ; kernel ; "
188
%
189
% Any extra ';' characters, at start, end or between kernel defintions are
190
% simply ignored.
191
%
192
% The special flags will expand a single kernel, into a list of rotated
193
% kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
194
% cyclic rotations, while a '>' will generate a list of 90-degree rotations.
195
% The '<' also exands using 90-degree rotates, but giving a 180-degree
196
% reflected kernel before the +/- 90-degree rotations, which can be important
197
% for Thinning operations.
198
%
199
167
% Note that 'name' kernels will start with an alphabetic character while the
200
168
% new kernel specification has a ':' character in its specification string.
201
169
% If neither is the case, it is assumed an old style of a simple list of
211
179
%
212
180
*/
213
181
214
/* This was separated so that it could be used as a separate
215
** array input handling function, such as for -color-matrix
216
*/
217
static KernelInfo *ParseKernelArray(const char *kernel_string)
182
MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
218
183
{
219
184
KernelInfo
220
185
*kernel;
222
187
char
223
188
token[MaxTextExtent];
224
189
190
register long
191
i;
192
225
193
const char
226
*p,
227
*end;
228
229
register ssize_t
230
i;
231
232
double
233
nan = sqrt((double)-1.0); /* Special Value : Not A Number */
194
*p;
234
195
235
196
MagickStatusType
236
197
flags;
238
199
GeometryInfo
239
200
args;
240
201
202
double
203
nan = sqrt((double)-1.0); /* Special Value : Not A Number */
204
205
assert(kernel_string != (const char *) NULL);
206
SetGeometryInfo(&args);
207
208
/* does it start with an alpha - Return a builtin kernel */
209
GetMagickToken(kernel_string,&p,token);
210
if ( isalpha((int)token[0]) )
211
{
212
long
213
type;
214
215
type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
216
if ( type < 0 || type == UserDefinedKernel )
217
return((KernelInfo *)NULL);
218
219
while (((isspace((int) ((unsigned char) *p)) != 0) ||
220
(*p == ',') || (*p == ':' )) && (*p != '\0'))
221
p++;
222
flags = ParseGeometry(p, &args);
223
224
/* special handling of missing values in input string */
225
switch( type ) {
226
case RectangleKernel:
227
if ( (flags & WidthValue) == 0 ) /* if no width then */
228
args.rho = args.sigma; /* then width = height */
229
if ( args.rho < 1.0 ) /* if width too small */
230
args.rho = 3; /* then width = 3 */
231
if ( args.sigma < 1.0 ) /* if height too small */
232
args.sigma = args.rho; /* then height = width */
233
if ( (flags & XValue) == 0 ) /* center offset if not defined */
234
args.xi = (double)(((long)args.rho-1)/2);
235
if ( (flags & YValue) == 0 )
236
args.psi = (double)(((long)args.sigma-1)/2);
237
break;
238
case SquareKernel:
239
case DiamondKernel:
240
case DiskKernel:
241
case PlusKernel:
242
if ( (flags & HeightValue) == 0 ) /* if no scale */
243
args.sigma = 1.0; /* then scale = 1.0 */
244
break;
245
default:
246
break;
247
}
248
249
return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
250
}
251
241
252
kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
242
253
if (kernel == (KernelInfo *)NULL)
243
254
return(kernel);
244
255
(void) ResetMagickMemory(kernel,0,sizeof(*kernel));
245
kernel->minimum = kernel->maximum = kernel->angle = 0.0;
246
kernel->negative_range = kernel->positive_range = 0.0;
247
256
kernel->type = UserDefinedKernel;
248
kernel->next = (KernelInfo *) NULL;
249
257
kernel->signature = MagickSignature;
250
258
251
/* find end of this specific kernel definition string */
252
end = strchr(kernel_string, ';');
253
if ( end == (char *) NULL )
254
end = strchr(kernel_string, '\0');
255
256
/* clear flags - for Expanding kernal lists thorugh rotations */
257
flags = NoValue;
258
259
259
/* Has a ':' in argument - New user kernel specification */
260
260
p = strchr(kernel_string, ':');
261
if ( p != (char *) NULL && p < end)
261
if ( p != (char *) NULL)
262
262
{
263
263
/* ParseGeometry() needs the geometry separated! -- Arrgghh */
264
264
memcpy(token, kernel_string, (size_t) (p-kernel_string));
265
265
token[p-kernel_string] = '\0';
266
SetGeometryInfo(&args);
267
266
flags = ParseGeometry(token, &args);
268
267
269
268
/* Size handling and checks of geometry settings */
273
272
args.rho = 1.0; /* then width = 1 */
274
273
if ( args.sigma < 1.0 ) /* if height too small */
275
274
args.sigma = args.rho; /* then height = width */
276
kernel->width = (size_t)args.rho;
277
kernel->height = (size_t)args.sigma;
275
kernel->width = (unsigned long)args.rho;
276
kernel->height = (unsigned long)args.sigma;
278
277
279
278
/* Offset Handling and Checks */
280
279
if ( args.xi < 0.0 || args.psi < 0.0 )
281
280
return(DestroyKernelInfo(kernel));
282
kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
283
: (ssize_t) (kernel->width-1)/2;
284
kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
285
: (ssize_t) (kernel->height-1)/2;
286
if ( kernel->x >= (ssize_t) kernel->width ||
287
kernel->y >= (ssize_t) kernel->height )
281
kernel->x = ((flags & XValue)!=0) ? (long)args.xi
282
: (long) (kernel->width-1)/2;
283
kernel->y = ((flags & YValue)!=0) ? (long)args.psi
284
: (long) (kernel->height-1)/2;
285
if ( kernel->x >= (long) kernel->width ||
286
kernel->y >= (long) kernel->height )
288
287
return(DestroyKernelInfo(kernel));
289
288
290
289
p++; /* advance beyond the ':' */
291
290
}
292
291
else
293
{ /* ELSE - Old old specification, forming odd-square kernel */
292
{ /* ELSE - Old old kernel specification, forming odd-square kernel */
294
293
/* count up number of values given */
295
294
p=(const char *) kernel_string;
296
295
while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
297
296
p++; /* ignore "'" chars for convolve filter usage - Cristy */
298
for (i=0; p < end; i++)
297
for (i=0; *p != '\0'; i++)
299
298
{
300
299
GetMagickToken(p,&p,token);
301
300
if (*token == ',')
302
301
GetMagickToken(p,&p,token);
303
302
}
304
303
/* set the size of the kernel - old sized square */
305
kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
306
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
304
kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
305
kernel->x = kernel->y = (long) (kernel->width-1)/2;
307
306
p=(const char *) kernel_string;
308
307
while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
309
308
p++; /* ignore "'" chars for convolve filter usage - Cristy */
318
317
kernel->minimum = +MagickHuge;
319
318
kernel->maximum = -MagickHuge;
320
319
kernel->negative_range = kernel->positive_range = 0.0;
321
322
for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
320
for (i=0; (i < (long) (kernel->width*kernel->height)) && (*p != '\0'); i++)
323
321
{
324
322
GetMagickToken(p,&p,token);
325
323
if (*token == ',')
326
324
GetMagickToken(p,&p,token);
327
325
if ( LocaleCompare("nan",token) == 0
328
|| LocaleCompare("-",token) == 0 ) {
326
|| LocaleCompare("-",token) == 0 ) {
329
327
kernel->values[i] = nan; /* do not include this value in kernel */
330
328
}
331
329
else {
337
335
Maximize(kernel->maximum, kernel->values[i]);
338
336
}
339
337
}
340
341
/* sanity check -- no more values in kernel definition */
342
GetMagickToken(p,&p,token);
343
if ( *token != '\0' && *token != ';' && *token != '\'' )
338
/* check that we recieved at least one real (non-nan) value! */
339
if ( kernel->minimum == MagickHuge )
344
340
return(DestroyKernelInfo(kernel));
345
341
346
#if 0
347
/* this was the old method of handling a incomplete kernel */
348
if ( i < (ssize_t) (kernel->width*kernel->height) ) {
342
/* This should not be needed for a fully defined kernel
343
* Perhaps an error should be reported instead!
344
* Kept for backward compatibility.
345
*/
346
if ( i < (long) (kernel->width*kernel->height) ) {
349
347
Minimize(kernel->minimum, kernel->values[i]);
350
348
Maximize(kernel->maximum, kernel->values[i]);
351
for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
349
for ( ; i < (long) (kernel->width*kernel->height); i++)
352
350
kernel->values[i]=0.0;
353
351
}
354
#else
355
/* Number of values for kernel was not enough - Report Error */
356
if ( i < (ssize_t) (kernel->width*kernel->height) )
357
return(DestroyKernelInfo(kernel));
358
#endif
359
360
/* check that we recieved at least one real (non-nan) value! */
361
if ( kernel->minimum == MagickHuge )
362
return(DestroyKernelInfo(kernel));
363
364
if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
365
ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
366
else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
367
ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
368
else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
369
ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
370
371
return(kernel);
372
}
373
374
static KernelInfo *ParseKernelName(const char *kernel_string)
375
{
376
KernelInfo
377
*kernel;
378
379
char
380
token[MaxTextExtent];
381
382
ssize_t
383
type;
384
385
const char
386
*p,
387
*end;
388
389
MagickStatusType
390
flags;
391
392
GeometryInfo
393
args;
394
395
/* Parse special 'named' kernel */
396
GetMagickToken(kernel_string,&p,token);
397
type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
398
if ( type < 0 || type == UserDefinedKernel )
399
return((KernelInfo *)NULL); /* not a valid named kernel */
400
401
while (((isspace((int) ((unsigned char) *p)) != 0) ||
402
(*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
403
p++;
404
405
end = strchr(p, ';'); /* end of this kernel defintion */
406
if ( end == (char *) NULL )
407
end = strchr(p, '\0');
408
409
/* ParseGeometry() needs the geometry separated! -- Arrgghh */
410
memcpy(token, p, (size_t) (end-p));
411
token[end-p] = '\0';
412
SetGeometryInfo(&args);
413
flags = ParseGeometry(token, &args);
414
415
#if 0
416
/* For Debugging Geometry Input */
417
fprintf(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
418
flags, args.rho, args.sigma, args.xi, args.psi );
419
#endif
420
421
/* special handling of missing values in input string */
422
switch( type ) {
423
case RectangleKernel:
424
if ( (flags & WidthValue) == 0 ) /* if no width then */
425
args.rho = args.sigma; /* then width = height */
426
if ( args.rho < 1.0 ) /* if width too small */
427
args.rho = 3; /* then width = 3 */
428
if ( args.sigma < 1.0 ) /* if height too small */
429
args.sigma = args.rho; /* then height = width */
430
if ( (flags & XValue) == 0 ) /* center offset if not defined */
431
args.xi = (double)(((ssize_t)args.rho-1)/2);
432
if ( (flags & YValue) == 0 )
433
args.psi = (double)(((ssize_t)args.sigma-1)/2);
434
break;
435
case SquareKernel:
436
case DiamondKernel:
437
case DiskKernel:
438
case PlusKernel:
439
case CrossKernel:
440
/* If no scale given (a 0 scale is valid! - set it to 1.0 */
441
if ( (flags & HeightValue) == 0 )
442
args.sigma = 1.0;
443
break;
444
case RingKernel:
445
if ( (flags & XValue) == 0 )
446
args.xi = 1.0;
447
break;
448
case ChebyshevKernel:
449
case ManhattanKernel:
450
case EuclideanKernel:
451
if ( (flags & HeightValue) == 0 ) /* no distance scale */
452
args.sigma = 100.0; /* default distance scaling */
453
else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
454
args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
455
else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
456
args.sigma *= QuantumRange/100.0; /* percentage of color range */
457
break;
458
default:
459
break;
460
}
461
462
kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
463
464
/* global expand to rotated kernel list - only for single kernels */
465
if ( kernel->next == (KernelInfo *) NULL ) {
466
if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
467
ExpandRotateKernelInfo(kernel, 45.0);
468
else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
469
ExpandRotateKernelInfo(kernel, 90.0);
470
else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
471
ExpandMirrorKernelInfo(kernel);
472
}
473
474
return(kernel);
475
}
476
477
MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
478
{
479
480
KernelInfo
481
*kernel,
482
*new_kernel;
483
484
char
485
token[MaxTextExtent];
486
487
const char
488
*p;
489
490
size_t
491
kernel_number;
492
493
p = kernel_string;
494
kernel = NULL;
495
kernel_number = 0;
496
497
while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
498
499
/* ignore extra or multiple ';' kernel seperators */
500
if ( *token != ';' ) {
501
502
/* tokens starting with alpha is a Named kernel */
503
if (isalpha((int) *token) != 0)
504
new_kernel = ParseKernelName(p);
505
else /* otherwise a user defined kernel array */
506
new_kernel = ParseKernelArray(p);
507
508
/* Error handling -- this is not proper error handling! */
509
if ( new_kernel == (KernelInfo *) NULL ) {
510
fprintf(stderr, "Failed to parse kernel number #%.20g\n",(double)
511
kernel_number);
512
if ( kernel != (KernelInfo *) NULL )
513
kernel=DestroyKernelInfo(kernel);
514
return((KernelInfo *) NULL);
515
}
516
517
/* initialise or append the kernel list */
518
if ( kernel == (KernelInfo *) NULL )
519
kernel = new_kernel;
520
else
521
LastKernelInfo(kernel)->next = new_kernel;
522
}
523
524
/* look for the next kernel in list */
525
p = strchr(p, ';');
526
if ( p == (char *) NULL )
527
break;
528
p++;
529
530
}
531
return(kernel);
532
}
533
352
353
return(kernel);
354
}
534
355
535
356
/*
536
357
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
564
385
%
565
386
% Convolution Kernels
566
387
%
567
% Unity
568
% the No-Op kernel, also requivelent to Gaussian of sigma zero.
569
% Basically a 3x3 kernel of a 1 surrounded by zeros.
570
%
571
% Gaussian:{radius},{sigma}
572
% Generate a two-dimentional gaussian kernel, as used by -gaussian.
573
% The sigma for the curve is required. The resulting kernel is
574
% normalized,
575
%
576
% If 'sigma' is zero, you get a single pixel on a field of zeros.
388
% Gaussian "{radius},{sigma}"
389
% Generate a two-dimentional gaussian kernel, as used by -gaussian
390
% A sigma is required, (with the 'x'), due to historical reasons.
577
391
%
578
392
% NOTE: that the 'radius' is optional, but if provided can limit (clip)
579
393
% the final size of the resulting kernel to a square 2*radius+1 in size.
582
396
% radius will be determined so as to produce the best minimal error
583
397
% result, which is usally much larger than is normally needed.
584
398
%
585
% LoG:{radius},{sigma}
586
% "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
587
% The supposed ideal edge detection, zero-summing kernel.
588
%
589
% An alturnative to this kernel is to use a "DoG" with a sigma ratio of
590
% approx 1.6 (according to wikipedia).
591
%
592
% DoG:{radius},{sigma1},{sigma2}
593
% "Difference of Gaussians" Kernel.
594
% As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
595
% from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
596
% The result is a zero-summing kernel.
597
%
598
% Blur:{radius},{sigma}[,{angle}]
599
% Generates a 1 dimensional or linear gaussian blur, at the angle given
600
% (current restricted to orthogonal angles). If a 'radius' is given the
601
% kernel is clipped to a width of 2*radius+1. Kernel can be rotated
602
% by a 90 degree angle.
603
%
604
% If 'sigma' is zero, you get a single pixel on a field of zeros.
605
%
606
% Note that two convolutions with two "Blur" kernels perpendicular to
607
% each other, is equivelent to a far larger "Gaussian" kernel with the
608
% same sigma value, However it is much faster to apply. This is how the
609
% "-blur" operator actually works.
610
%
611
% Comet:{width},{sigma},{angle}
612
% Blur in one direction only, much like how a bright object leaves
399
% Blur "{radius},{sigma},{angle}"
400
% As per Gaussian, but generates a 1 dimensional or linear gaussian
401
% blur, at the angle given (current restricted to orthogonal angles).
402
% If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
403
%
404
% NOTE that two such blurs perpendicular to each other is equivelent to
405
% -blur and the previous gaussian, but is often 10 or more times faster.
406
%
407
% Comet "{width},{sigma},{angle}"
408
% Blur in one direction only, mush like how a bright object leaves
613
409
% a comet like trail. The Kernel is actually half a gaussian curve,
614
% Adding two such blurs in opposite directions produces a Blur Kernel.
615
% Angle can be rotated in multiples of 90 degrees.
410
% Adding two such blurs in oppiste directions produces a Linear Blur.
616
411
%
617
% Note that the first argument is the width of the kernel and not the
412
% NOTE: that the first argument is the width of the kernel and not the
618
413
% radius of the kernel.
619
414
%
620
415
% # Still to be implemented...
621
416
% #
417
% # Sharpen "{radius},{sigma}
418
% # Negated Gaussian (center zeroed and re-normalized),
419
% # with a 2 unit positive peak. -- Check On line documentation
420
% #
421
% # Laplacian "{radius},{sigma}"
422
% # Laplacian (a mexican hat like) Function
423
% #
424
% # LOG "{radius},{sigma1},{sigma2}
425
% # Laplacian of Gaussian
426
% #
427
% # DOG "{radius},{sigma1},{sigma2}
428
% # Difference of two Gaussians
429
% #
622
430
% # Filter2D
623
431
% # Filter1D
624
432
% # Set kernel values using a resize filter, and given scale (sigma)
625
433
% # Cylindrical or Linear. Is this posible with an image?
626
434
% #
627
435
%
628
% Named Constant Convolution Kernels
629
%
630
% All these are unscaled, zero-summing kernels by default. As such for
631
% non-HDRI version of ImageMagick some form of normalization, user scaling,
632
% and biasing the results is recommended, to prevent the resulting image
633
% being 'clipped'.
634
%
635
% The 3x3 kernels (most of these) can be circularly rotated in multiples of
636
% 45 degrees to generate the 8 angled varients of each of the kernels.
637
%
638
% Laplacian:{type}
639
% Discrete Lapacian Kernels, (without normalization)
640
% Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
641
% Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
642
% Type 2 : 3x3 with center:4 edge:1 corner:-2
643
% Type 3 : 3x3 with center:4 edge:-2 corner:1
644
% Type 5 : 5x5 laplacian
645
% Type 7 : 7x7 laplacian
646
% Type 15 : 5x5 LoG (sigma approx 1.4)
647
% Type 19 : 9x9 LoG (sigma approx 1.4)
648
%
649
% Sobel:{angle}
650
% Sobel 'Edge' convolution kernel (3x3)
651
% | -1, 0, 1 |
652
% | -2, 0,-2 |
653
% | -1, 0, 1 |
654
%
655
% Sobel:{type},{angle}
656
% Type 0: default un-nomalized version shown above.
657
%
658
% Type 1: As default but pre-normalized
659
% | 1, 0, -1 |
660
% | 2, 0, -2 | / 4
661
% | 1, 0, -1 |
662
%
663
% Type 2: Diagonal version with same normalization as 1
664
% | 1, 0, -1 |
665
% | 2, 0, -2 | / 4
666
% | 1, 0, -1 |
667
%
668
% Roberts:{angle}
669
% Roberts convolution kernel (3x3)
670
% | 0, 0, 0 |
671
% | -1, 1, 0 |
672
% | 0, 0, 0 |
673
%
674
% Prewitt:{angle}
675
% Prewitt Edge convolution kernel (3x3)
676
% | -1, 0, 1 |
677
% | -1, 0, 1 |
678
% | -1, 0, 1 |
679
%
680
% Compass:{angle}
681
% Prewitt's "Compass" convolution kernel (3x3)
682
% | -1, 1, 1 |
683
% | -1,-2, 1 |
684
% | -1, 1, 1 |
685
%
686
% Kirsch:{angle}
687
% Kirsch's "Compass" convolution kernel (3x3)
688
% | -3,-3, 5 |
689
% | -3, 0, 5 |
690
% | -3,-3, 5 |
691
%
692
% FreiChen:{angle}
693
% Frei-Chen Edge Detector is based on a kernel that is similar to
694
% the Sobel Kernel, but is designed to be isotropic. That is it takes
695
% into account the distance of the diagonal in the kernel.
696
%
697
% | 1, 0, -1 |
698
% | sqrt(2), 0, -sqrt(2) |
699
% | 1, 0, -1 |
700
%
701
% FreiChen:{type},{angle}
702
%
703
% Frei-Chen Pre-weighted kernels...
704
%
705
% Type 0: default un-nomalized version shown above.
706
%
707
% Type 1: Orthogonal Kernel (same as type 11 below)
708
% | 1, 0, -1 |
709
% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
710
% | 1, 0, -1 |
711
%
712
% Type 2: Diagonal form of Kernel...
713
% | 1, sqrt(2), 0 |
714
% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
715
% | 0, -sqrt(2) -1 |
716
%
717
% However this kernel is als at the heart of the FreiChen Edge Detection
718
% Process which uses a set of 9 specially weighted kernel. These 9
719
% kernels not be normalized, but directly applied to the image. The
720
% results is then added together, to produce the intensity of an edge in
721
% a specific direction. The square root of the pixel value can then be
722
% taken as the cosine of the edge, and at least 2 such runs at 90 degrees
723
% from each other, both the direction and the strength of the edge can be
724
% determined.
725
%
726
% Type 10: All 9 of the following pre-weighted kernels...
727
%
728
% Type 11: | 1, 0, -1 |
729
% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
730
% | 1, 0, -1 |
731
%
732
% Type 12: | 1, sqrt(2), 1 |
733
% | 0, 0, 0 | / 2*sqrt(2)
734
% | 1, sqrt(2), 1 |
735
%
736
% Type 13: | sqrt(2), -1, 0 |
737
% | -1, 0, 1 | / 2*sqrt(2)
738
% | 0, 1, -sqrt(2) |
739
%
740
% Type 14: | 0, 1, -sqrt(2) |
741
% | -1, 0, 1 | / 2*sqrt(2)
742
% | sqrt(2), -1, 0 |
743
%
744
% Type 15: | 0, -1, 0 |
745
% | 1, 0, 1 | / 2
746
% | 0, -1, 0 |
747
%
748
% Type 16: | 1, 0, -1 |
749
% | 0, 0, 0 | / 2
750
% | -1, 0, 1 |
751
%
752
% Type 17: | 1, -2, 1 |
753
% | -2, 4, -2 | / 6
754
% | -1, -2, 1 |
755
%
756
% Type 18: | -2, 1, -2 |
757
% | 1, 4, 1 | / 6
758
% | -2, 1, -2 |
759
%
760
% Type 19: | 1, 1, 1 |
761
% | 1, 1, 1 | / 3
762
% | 1, 1, 1 |
763
%
764
% The first 4 are for edge detection, the next 4 are for line detection
765
% and the last is to add a average component to the results.
766
%
767
% Using a special type of '-1' will return all 9 pre-weighted kernels
768
% as a multi-kernel list, so that you can use them directly (without
769
% normalization) with the special "-set option:morphology:compose Plus"
770
% setting to apply the full FreiChen Edge Detection Technique.
771
%
772
% If 'type' is large it will be taken to be an actual rotation angle for
773
% the default FreiChen (type 0) kernel. As such FreiChen:45 will look
774
% like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
775
%
776
% WARNING: The above was layed out as per
777
% http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
778
% But rotated 90 degrees so direction is from left rather than the top.
779
% I have yet to find any secondary confirmation of the above. The only
780
% other source found was actual source code at
781
% http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
782
% Neigher paper defineds the kernels in a way that looks locical or
783
% correct when taken as a whole.
784
%
785
436
% Boolean Kernels
786
437
%
787
% Diamond:[{radius}[,{scale}]]
788
% Generate a diamond shaped kernel with given radius to the points.
438
% Rectangle "{geometry}"
439
% Simply generate a rectangle of 1's with the size given. You can also
440
% specify the location of the 'control point', otherwise the closest
441
% pixel to the center of the rectangle is selected.
442
%
443
% Properly centered and odd sized rectangles work the best.
444
%
445
% Diamond "[{radius}[,{scale}]]"
446
% Generate a diamond shaped kernal with given radius to the points.
789
447
% Kernel size will again be radius*2+1 square and defaults to radius 1,
790
448
% generating a 3x3 kernel that is slightly larger than a square.
791
449
%
792
% Square:[{radius}[,{scale}]]
450
% Square "[{radius}[,{scale}]]"
793
451
% Generate a square shaped kernel of size radius*2+1, and defaulting
794
452
% to a 3x3 (radius 1).
795
453
%
796
% Note that using a larger radius for the "Square" or the "Diamond" is
797
% also equivelent to iterating the basic morphological method that many
798
% times. However iterating with the smaller radius is actually faster
799
% than using a larger kernel radius.
800
%
801
% Rectangle:{geometry}
802
% Simply generate a rectangle of 1's with the size given. You can also
803
% specify the location of the 'control point', otherwise the closest
804
% pixel to the center of the rectangle is selected.
805
%
806
% Properly centered and odd sized rectangles work the best.
807
%
808
% Disk:[{radius}[,{scale}]]
454
% Note that using a larger radius for the "Square" or the "Diamond"
455
% is also equivelent to iterating the basic morphological method
456
% that many times. However However iterating with the smaller radius 1
457
% default is actually faster than using a larger kernel radius.
458
%
459
% Disk "[{radius}[,{scale}]]
809
460
% Generate a binary disk of the radius given, radius may be a float.
810
461
% Kernel size will be ceil(radius)*2+1 square.
811
462
% NOTE: Here are some disk shapes of specific interest
812
% "Disk:1" => "diamond" or "cross:1"
813
% "Disk:1.5" => "square"
814
% "Disk:2" => "diamond:2"
815
% "Disk:2.5" => a general disk shape of radius 2
816
% "Disk:2.9" => "square:2"
817
% "Disk:3.5" => default - octagonal/disk shape of radius 3
818
% "Disk:4.2" => roughly octagonal shape of radius 4
819
% "Disk:4.3" => a general disk shape of radius 4
463
% "disk:1" => "diamond" or "cross:1"
464
% "disk:1.5" => "square"
465
% "disk:2" => "diamond:2"
466
% "disk:2.5" => a general disk shape of radius 2
467
% "disk:2.9" => "square:2"
468
% "disk:3.5" => default - octagonal/disk shape of radius 3
469
% "disk:4.2" => roughly octagonal shape of radius 4
470
% "disk:4.3" => a general disk shape of radius 4
820
471
% After this all the kernel shape becomes more and more circular.
821
472
%
822
473
% Because a "disk" is more circular when using a larger radius, using a
823
474
% larger radius is preferred over iterating the morphological operation.
824
475
%
825
% Symbol Dilation Kernels
826
%
827
% These kernel is not a good general morphological kernel, but is used
828
% more for highlighting and marking any single pixels in an image using,
829
% a "Dilate" method as appropriate.
830
%
831
% For the same reasons iterating these kernels does not produce the
832
% same result as using a larger radius for the symbol.
833
%
834
% Plus:[{radius}[,{scale}]]
835
% Cross:[{radius}[,{scale}]]
836
% Generate a kernel in the shape of a 'plus' or a 'cross' with
837
% a each arm the length of the given radius (default 2).
476
% Plus "[{radius}[,{scale}]]"
477
% Generate a kernel in the shape of a 'plus' sign. The length of each
478
% arm is also the radius, which defaults to 2.
479
%
480
% This kernel is not a good general morphological kernel, but is used
481
% more for highlighting and marking any single pixels in an image using,
482
% a "Dilate" or "Erode" method as appropriate.
838
483
%
839
484
% NOTE: "plus:1" is equivelent to a "Diamond" kernel.
840
485
%
841
% Ring:{radius1},{radius2}[,{scale}]
842
% A ring of the values given that falls between the two radii.
843
% Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
844
% This is the 'edge' pixels of the default "Disk" kernel,
845
% More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
846
%
847
% Hit and Miss Kernels
848
%
849
% Peak:radius1,radius2
850
% Find any peak larger than the pixels the fall between the two radii.
851
% The default ring of pixels is as per "Ring".
852
% Edges
853
% Find flat orthogonal edges of a binary shape
854
% Corners
855
% Find 90 degree corners of a binary shape
856
% LineEnds:type
857
% Find end points of lines (for pruning a skeletion)
858
% Two types of lines ends (default to both) can be searched for
859
% Type 0: All line ends
860
% Type 1: single kernel for 4-conneected line ends
861
% Type 2: single kernel for simple line ends
862
% LineJunctions
863
% Find three line junctions (within a skeletion)
864
% Type 0: all line junctions
865
% Type 1: Y Junction kernel
866
% Type 2: Diagonal T Junction kernel
867
% Type 3: Orthogonal T Junction kernel
868
% Type 4: Diagonal X Junction kernel
869
% Type 5: Orthogonal + Junction kernel
870
% Ridges:type
871
% Find single pixel ridges or thin lines
872
% Type 1: Fine single pixel thick lines and ridges
873
% Type 2: Find two pixel thick lines and ridges
874
% ConvexHull
875
% Octagonal thicken kernel, to generate convex hulls of 45 degrees
876
% Skeleton:type
877
% Traditional skeleton generating kernels.
878
% Type 1: Tradional Skeleton kernel (4 connected skeleton)
879
% Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
880
% Type 3: Experimental Variation to try to present left-right symmetry
881
% Type 4: Experimental Variation to preserve left-right symmetry
486
% Note that unlike other kernels iterating a plus does not produce the
487
% same result as using a larger radius for the cross.
882
488
%
883
489
% Distance Measuring Kernels
884
490
%
885
% Different types of distance measuring methods, which are used with the
886
% a 'Distance' morphology method for generating a gradient based on
887
% distance from an edge of a binary shape, though there is a technique
888
% for handling a anti-aliased shape.
889
%
890
% See the 'Distance' Morphological Method, for information of how it is
891
% applied.
892
%
893
% Chebyshev:[{radius}][x{scale}[%!]]
491
% Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1)
492
% Manhatten "[{radius}][x{scale}]" square grid distance (default r=1)
493
% Euclidean "[{radius}][x{scale}]" direct distance (default r=1)
494
%
495
% Different types of distance measuring methods, which are used with the
496
% a 'Distance' morphology method for generating a gradient based on
497
% distance from an edge of a binary shape, though there is a technique
498
% for handling a anti-aliased shape.
499
%
894
500
% Chebyshev Distance (also known as Tchebychev Distance) is a value of
895
501
% one to any neighbour, orthogonal or diagonal. One why of thinking of
896
502
% it is the number of squares a 'King' or 'Queen' in chess needs to
898
504
% 'square' like distance function, but one where diagonals are closer
899
505
% than expected.
900
506
%
901
% Manhattan:[{radius}][x{scale}[%!]]
902
% Manhattan Distance (also known as Rectilinear Distance, or the Taxi
507
% Manhatten Distance (also known as Rectilinear Distance, or the Taxi
903
508
% Cab metric), is the distance needed when you can only travel in
904
509
% orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
905
510
% in chess would travel. It results in a diamond like distances, where
906
511
% diagonals are further than expected.
907
512
%
908
% Euclidean:[{radius}][x{scale}[%!]]
909
513
% Euclidean Distance is the 'direct' or 'as the crow flys distance.
910
514
% However by default the kernel size only has a radius of 1, which
911
515
% limits the distance to 'Knight' like moves, with only orthogonal and
920
524
% To allow the use of fractional distances that you get with diagonals
921
525
% the actual distance is scaled by a fixed value which the user can
922
526
% provide. This is not actually nessary for either ""Chebyshev" or
923
% "Manhattan" distance kernels, but is done for all three distance
527
% "Manhatten" distance kernels, but is done for all three distance
924
528
% kernels. If no scale is provided it is set to a value of 100,
925
529
% allowing for a maximum distance measurement of 655 pixels using a Q16
926
530
% version of IM, from any edge. However for small images this can
927
531
% result in quite a dark gradient.
928
532
%
533
% See the 'Distance' Morphological Method, for information of how it is
534
% applied.
535
%
536
% # Hit-n-Miss Kernel-Lists -- Still to be implemented
537
% #
538
% # specifically for Pruning, Thinning, Thickening
539
% #
929
540
*/
930
541
931
542
MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
934
545
KernelInfo
935
546
*kernel;
936
547
937
register ssize_t
548
register long
938
549
i;
939
550
940
register ssize_t
551
register long
941
552
u,
942
553
v;
943
554
944
555
double
945
556
nan = sqrt((double)-1.0); /* Special Value : Not A Number */
946
557
947
/* Generate a new empty kernel if needed */
948
kernel=(KernelInfo *) NULL;
949
switch(type) {
950
case UndefinedKernel: /* These should not call this function */
951
case UserDefinedKernel:
952
break;
953
case UnityKernel: /* Named Descrete Convolution Kernels */
954
case LaplacianKernel:
955
case SobelKernel:
956
case RobertsKernel:
957
case PrewittKernel:
958
case CompassKernel:
959
case KirschKernel:
960
case FreiChenKernel:
961
case EdgesKernel: /* Hit and Miss kernels */
962
case CornersKernel:
963
case ThinDiagonalsKernel:
964
case LineEndsKernel:
965
case LineJunctionsKernel:
966
case RidgesKernel:
967
case ConvexHullKernel:
968
case SkeletonKernel:
969
break; /* A pre-generated kernel is not needed */
970
#if 0
971
/* set to 1 to do a compile-time check that we haven't missed anything */
972
case GaussianKernel:
973
case DoGKernel:
974
case LoGKernel:
975
case BlurKernel:
976
case CometKernel:
977
case DiamondKernel:
978
case SquareKernel:
979
case RectangleKernel:
980
case DiskKernel:
981
case PlusKernel:
982
case CrossKernel:
983
case RingKernel:
984
case PeaksKernel:
985
case ChebyshevKernel:
986
case ManhattanKernel:
987
case EuclideanKernel:
988
#else
989
default:
990
#endif
991
/* Generate the base Kernel Structure */
992
kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
993
if (kernel == (KernelInfo *) NULL)
994
return(kernel);
995
(void) ResetMagickMemory(kernel,0,sizeof(*kernel));
996
kernel->minimum = kernel->maximum = kernel->angle = 0.0;
997
kernel->negative_range = kernel->positive_range = 0.0;
998
kernel->type = type;
999
kernel->next = (KernelInfo *) NULL;
1000
kernel->signature = MagickSignature;
1001
break;
1002
}
558
kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
559
if (kernel == (KernelInfo *) NULL)
560
return(kernel);
561
(void) ResetMagickMemory(kernel,0,sizeof(*kernel));
562
kernel->minimum = kernel->maximum = 0.0;
563
kernel->negative_range = kernel->positive_range = 0.0;
564
kernel->type = type;
565
kernel->signature = MagickSignature;
1003
566
1004
567
switch(type) {
1005
568
/* Convolution Kernels */
1006
569
case GaussianKernel:
1007
case DoGKernel:
1008
case LoGKernel:
1009
570
{ double
1010
sigma = fabs(args->sigma),
1011
sigma2 = fabs(args->xi),
1012
A, B, R;
1013
1014
if ( args->rho >= 1.0 )
1015
kernel->width = (size_t)args->rho*2+1;
1016
else if ( (type != DoGKernel) || (sigma >= sigma2) )
1017
kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1018
else
1019
kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1020
kernel->height = kernel->width;
1021
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
571
sigma = fabs(args->sigma);
572
573
sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
574
575
kernel->width = kernel->height =
576
GetOptimalKernelWidth2D(args->rho,sigma);
577
kernel->x = kernel->y = (long) (kernel->width-1)/2;
578
kernel->negative_range = kernel->positive_range = 0.0;
1022
579
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1023
580
kernel->height*sizeof(double));
1024
581
if (kernel->values == (double *) NULL)
1025
582
return(DestroyKernelInfo(kernel));
1026
583
1027
/* WARNING: The following generates a 'sampled gaussian' kernel.
1028
* What we really want is a 'discrete gaussian' kernel.
1029
*
1030
* How to do this is currently not known, but appears to be
1031
* basied on the Error Function 'erf()' (intergral of a gaussian)
1032
*/
1033
1034
if ( type == GaussianKernel || type == DoGKernel )
1035
{ /* Calculate a Gaussian, OR positive half of a DoG */
1036
if ( sigma > MagickEpsilon )
1037
{ A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1038
B = 1.0/(Magick2PI*sigma*sigma);
1039
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1040
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1041
kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1042
}
1043
else /* limiting case - a unity (normalized Dirac) kernel */
1044
{ (void) ResetMagickMemory(kernel->values,0, (size_t)
1045
kernel->width*kernel->height*sizeof(double));
1046
kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1047
}
1048
}
1049
1050
if ( type == DoGKernel )
1051
{ /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1052
if ( sigma2 > MagickEpsilon )
1053
{ sigma = sigma2; /* simplify loop expressions */
1054
A = 1.0/(2.0*sigma*sigma);
1055
B = 1.0/(Magick2PI*sigma*sigma);
1056
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1057
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1058
kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1059
}
1060
else /* limiting case - a unity (normalized Dirac) kernel */
1061
kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1062
}
1063
1064
if ( type == LoGKernel )
1065
{ /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1066
if ( sigma > MagickEpsilon )
1067
{ A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1068
B = 1.0/(MagickPI*sigma*sigma*sigma*sigma);
1069
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1070
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1071
{ R = ((double)(u*u+v*v))*A;
1072
kernel->values[i] = (1-R)*exp(-R)*B;
1073
}
1074
}
1075
else /* special case - generate a unity kernel */
1076
{ (void) ResetMagickMemory(kernel->values,0, (size_t)
1077
kernel->width*kernel->height*sizeof(double));
1078
kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1079
}
1080
}
1081
1082
/* Note the above kernels may have been 'clipped' by a user defined
1083
** radius, producing a smaller (darker) kernel. Also for very small
1084
** sigma's (> 0.1) the central value becomes larger than one, and thus
1085
** producing a very bright kernel.
1086
**
1087
** Normalization will still be needed.
1088
*/
1089
1090
/* Normalize the 2D Gaussian Kernel
1091
**
1092
** NB: a CorrelateNormalize performs a normal Normalize if
1093
** there are no negative values.
1094
*/
1095
CalcKernelMetaData(kernel); /* the other kernel meta-data */
1096
ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
584
sigma = 2.0*sigma*sigma; /* simplify the expression */
585
for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
586
for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
587
kernel->positive_range += (
588
kernel->values[i] =
589
exp(-((double)(u*u+v*v))/sigma)
590
/* / (MagickPI*sigma) */ );
591
kernel->minimum = 0;
592
kernel->maximum = kernel->values[
593
kernel->y*kernel->width+kernel->x ];
594
595
ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1097
596
1098
597
break;
1099
598
}
1100
599
case BlurKernel:
1101
600
{ double
1102
sigma = fabs(args->sigma),
1103
alpha, beta;
1104
1105
if ( args->rho >= 1.0 )
1106
kernel->width = (size_t)args->rho*2+1;
1107
else
1108
kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
601
sigma = fabs(args->sigma);
602
603
sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
604
605
kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
606
kernel->x = (long) (kernel->width-1)/2;
1109
607
kernel->height = 1;
1110
kernel->x = (ssize_t) (kernel->width-1)/2;
1111
608
kernel->y = 0;
1112
609
kernel->negative_range = kernel->positive_range = 0.0;
1113
610
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1124
621
** As such while wierd it is prefered.
1125
622
**
1126
623
** I am told this method originally came from Photoshop.
1127
**
1128
** A properly normalized curve is generated (apart from edge clipping)
1129
** even though we later normalize the result (for edge clipping)
1130
** to allow the correct generation of a "Difference of Blurs".
1131
624
*/
1132
1133
/* initialize */
1134
v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
625
sigma *= KernelRank; /* simplify expanded curve */
626
v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1135
627
(void) ResetMagickMemory(kernel->values,0, (size_t)
1136
kernel->width*kernel->height*sizeof(double));
1137
/* Calculate a Positive 1D Gaussian */
1138
if ( sigma > MagickEpsilon )
1139
{ sigma *= KernelRank; /* simplify loop expressions */
1140
alpha = 1.0/(2.0*sigma*sigma);
1141
beta= 1.0/(MagickSQ2PI*sigma );
1142
for ( u=-v; u <= v; u++) {
1143
kernel->values[(u+v)/KernelRank] +=
1144
exp(-((double)(u*u))*alpha)*beta;
1145
}
1146
}
1147
else /* special case - generate a unity kernel */
1148
kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
628
kernel->width*sizeof(double));
629
for ( u=-v; u <= v; u++) {
630
kernel->values[(u+v)/KernelRank] +=
631
exp(-((double)(u*u))/(2.0*sigma*sigma))
632
/* / (MagickSQ2PI*sigma/KernelRank) */ ;
633
}
634
for (i=0; i < (long) kernel->width; i++)
635
kernel->positive_range += kernel->values[i];
1149
636
#else
1150
/* Direct calculation without curve averaging */
1151
1152
/* Calculate a Positive Gaussian */
1153
if ( sigma > MagickEpsilon )
1154
{ alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1155
beta = 1.0/(MagickSQ2PI*sigma);
1156
for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1157
kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1158
}
1159
else /* special case - generate a unity kernel */
1160
{ (void) ResetMagickMemory(kernel->values,0, (size_t)
1161
kernel->width*kernel->height*sizeof(double));
1162
kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1163
}
637
for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
638
kernel->positive_range += (
639
kernel->values[i] =
640
exp(-((double)(u*u))/(2.0*sigma*sigma))
641
/* / (MagickSQ2PI*sigma) */ );
1164
642
#endif
1165
/* Note the above kernel may have been 'clipped' by a user defined
643
kernel->minimum = 0;
644
kernel->maximum = kernel->values[ kernel->x ];
645
/* Note that neither methods above generate a normalized kernel,
646
** though it gets close. The kernel may be 'clipped' by a user defined
1166
647
** radius, producing a smaller (darker) kernel. Also for very small
1167
648
** sigma's (> 0.1) the central value becomes larger than one, and thus
1168
649
** producing a very bright kernel.
1169
**
1170
** Normalization will still be needed.
1171
650
*/
1172
651
1173
652
/* Normalize the 1D Gaussian Kernel
1174
653
**
1175
** NB: a CorrelateNormalize performs a normal Normalize if
1176
** there are no negative values.
654
** Because of this the divisor in the above kernel generator is
655
** not needed, so is not done above.
1177
656
*/
1178
CalcKernelMetaData(kernel); /* the other kernel meta-data */
1179
ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
657
ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1180
658
1181
/* rotate the 1D kernel by given angle */
1182
RotateKernelInfo(kernel, args->xi );
659
/* rotate the kernel by given angle */
660
RotateKernelInfo(kernel, args->xi);
1183
661
break;
1184
662
}
1185
663
case CometKernel:
1186
664
{ double
1187
sigma = fabs(args->sigma),
1188
A;
665
sigma = fabs(args->sigma);
666
667
sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
1189
668
1190
669
if ( args->rho < 1.0 )
1191
kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
670
kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1192
671
else
1193
kernel->width = (size_t)args->rho;
672
kernel->width = (unsigned long)args->rho;
1194
673
kernel->x = kernel->y = 0;
1195
674
kernel->height = 1;
1196
675
kernel->negative_range = kernel->positive_range = 0.0;
1199
678
if (kernel->values == (double *) NULL)
1200
679
return(DestroyKernelInfo(kernel));
1201
680
1202
/* A comet blur is half a 1D gaussian curve, so that the object is
681
/* A comet blur is half a gaussian curve, so that the object is
1203
682
** blurred in one direction only. This may not be quite the right
1204
** curve to use so may change in the future. The function must be
1205
** normalised after generation, which also resolves any clipping.
1206
**
1207
** As we are normalizing and not subtracting gaussians,
1208
** there is no need for a divisor in the gaussian formula
1209
**
1210
** It is less comples
683
** curve so may change in the future. The function must be normalised.
1211
684
*/
1212
if ( sigma > MagickEpsilon )
1213
{
1214
685
#if 1
1215
686
#define KernelRank 3
1216
v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1217
(void) ResetMagickMemory(kernel->values,0, (size_t)
1218
kernel->width*sizeof(double));
1219
sigma *= KernelRank; /* simplify the loop expression */
1220
A = 1.0/(2.0*sigma*sigma);
1221
/* B = 1.0/(MagickSQ2PI*sigma); */
1222
for ( u=0; u < v; u++) {
1223
kernel->values[u/KernelRank] +=
1224
exp(-((double)(u*u))*A);
1225
/* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1226
}
1227
for (i=0; i < (ssize_t) kernel->width; i++)
1228
kernel->positive_range += kernel->values[i];
687
sigma *= KernelRank; /* simplify expanded curve */
688
v = (long) kernel->width*KernelRank; /* start/end points to fit range */
689
(void) ResetMagickMemory(kernel->values,0, (size_t)
690
kernel->width*sizeof(double));
691
for ( u=0; u < v; u++) {
692
kernel->values[u/KernelRank] +=
693
exp(-((double)(u*u))/(2.0*sigma*sigma))
694
/* / (MagickSQ2PI*sigma/KernelRank) */ ;
695
}
696
for (i=0; i < (long) kernel->width; i++)
697
kernel->positive_range += kernel->values[i];
1229
698
#else
1230
A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1231
/* B = 1.0/(MagickSQ2PI*sigma); */
1232
for ( i=0; i < (ssize_t) kernel->width; i++)
1233
kernel->positive_range +=
1234
kernel->values[i] =
1235
exp(-((double)(i*i))*A);
1236
/* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
699
for ( i=0; i < (long) kernel->width; i++)
700
kernel->positive_range += (
701
kernel->values[i] =
702
exp(-((double)(i*i))/(2.0*sigma*sigma))
703
/* / (MagickSQ2PI*sigma) */ );
1237
704
#endif
1238
}
1239
else /* special case - generate a unity kernel */
1240
{ (void) ResetMagickMemory(kernel->values,0, (size_t)
1241
kernel->width*kernel->height*sizeof(double));
1242
kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1243
kernel->positive_range = 1.0;
1244
}
1245
1246
kernel->minimum = 0.0;
705
kernel->minimum = 0;
1247
706
kernel->maximum = kernel->values[0];
1248
kernel->negative_range = 0.0;
1249
707
1250
708
ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1251
709
RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1252
710
break;
1253
711
}
1254
1255
/* Convolution Kernels - Well Known Constants */
1256
case LaplacianKernel:
1257
{ switch ( (int) args->rho ) {
1258
case 0:
1259
default: /* laplacian square filter -- default */
1260
kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1261
break;
1262
case 1: /* laplacian diamond filter */
1263
kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1264
break;
1265
case 2:
1266
kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1267
break;
1268
case 3:
1269
kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1270
break;
1271
case 5: /* a 5x5 laplacian */
1272
kernel=ParseKernelArray(
1273
"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");
1274
break;
1275
case 7: /* a 7x7 laplacian */
1276
kernel=ParseKernelArray(
1277
"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" );
1278
break;
1279
case 15: /* a 5x5 LoG (sigma approx 1.4) */
1280
kernel=ParseKernelArray(
1281
"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");
1282
break;
1283
case 19: /* a 9x9 LoG (sigma approx 1.4) */
1284
/* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1285
kernel=ParseKernelArray(
1286
"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");
1287
break;
1288
}
1289
if (kernel == (KernelInfo *) NULL)
1290
return(kernel);
1291
kernel->type = type;
1292
break;
1293
}
1294
case SobelKernel:
1295
#if 0
1296
{ /* Sobel with optional 'sub-types' */
1297
switch ( (int) args->rho ) {
1298
default:
1299
case 0:
1300
kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1301
if (kernel == (KernelInfo *) NULL)
1302
return(kernel);
1303
kernel->type = type;
1304
break;
1305
case 1:
1306
kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1307
if (kernel == (KernelInfo *) NULL)
1308
return(kernel);
1309
kernel->type = type;
1310
ScaleKernelInfo(kernel, 0.25, NoValue);
1311
break;
1312
case 2:
1313
kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1314
if (kernel == (KernelInfo *) NULL)
1315
return(kernel);
1316
kernel->type = type;
1317
ScaleKernelInfo(kernel, 0.25, NoValue);
1318
break;
1319
}
1320
if ( fabs(args->sigma) > MagickEpsilon )
1321
/* Rotate by correctly supplied 'angle' */
1322
RotateKernelInfo(kernel, args->sigma);
1323
else if ( args->rho > 30.0 || args->rho < -30.0 )
1324
/* Rotate by out of bounds 'type' */
1325
RotateKernelInfo(kernel, args->rho);
1326
break;
1327
}
1328
#else
1329
{ /* Simple Sobel Kernel */
1330
kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1331
if (kernel == (KernelInfo *) NULL)
1332
return(kernel);
1333
kernel->type = type;
1334
RotateKernelInfo(kernel, args->rho);
1335
break;
1336
}
1337
#endif
1338
case RobertsKernel:
1339
{
1340
kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1341
if (kernel == (KernelInfo *) NULL)
1342
return(kernel);
1343
kernel->type = type;
1344
RotateKernelInfo(kernel, args->rho);
1345
break;
1346
}
1347
case PrewittKernel:
1348
{
1349
kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1350
if (kernel == (KernelInfo *) NULL)
1351
return(kernel);
1352
kernel->type = type;
1353
RotateKernelInfo(kernel, args->rho);
1354
break;
1355
}
1356
case CompassKernel:
1357
{
1358
kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1359
if (kernel == (KernelInfo *) NULL)
1360
return(kernel);
1361
kernel->type = type;
1362
RotateKernelInfo(kernel, args->rho);
1363
break;
1364
}
1365
case KirschKernel:
1366
{
1367
kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1368
if (kernel == (KernelInfo *) NULL)
1369
return(kernel);
1370
kernel->type = type;
1371
RotateKernelInfo(kernel, args->rho);
1372
break;
1373
}
1374
case FreiChenKernel:
1375
/* Direction is set to be left to right positive */
1376
/* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1377
/* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1378
{ switch ( (int) args->rho ) {
1379
default:
1380
case 0:
1381
kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1382
if (kernel == (KernelInfo *) NULL)
1383
return(kernel);
1384
kernel->type = type;
1385
kernel->values[3] = +MagickSQ2;
1386
kernel->values[5] = -MagickSQ2;
1387
CalcKernelMetaData(kernel); /* recalculate meta-data */
1388
break;
1389
case 2:
1390
kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1391
if (kernel == (KernelInfo *) NULL)
1392
return(kernel);
1393
kernel->type = type;
1394
kernel->values[1] = kernel->values[3] = +MagickSQ2;
1395
kernel->values[5] = kernel->values[7] = -MagickSQ2;
1396
CalcKernelMetaData(kernel); /* recalculate meta-data */
1397
ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1398
break;
1399
case 10:
1400
kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1401
if (kernel == (KernelInfo *) NULL)
1402
return(kernel);
1403
break;
1404
case 1:
1405
case 11:
1406
kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1407
if (kernel == (KernelInfo *) NULL)
1408
return(kernel);
1409
kernel->type = type;
1410
kernel->values[3] = +MagickSQ2;
1411
kernel->values[5] = -MagickSQ2;
1412
CalcKernelMetaData(kernel); /* recalculate meta-data */
1413
ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1414
break;
1415
case 12:
1416
kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1417
if (kernel == (KernelInfo *) NULL)
1418
return(kernel);
1419
kernel->type = type;
1420
kernel->values[1] = +MagickSQ2;
1421
kernel->values[7] = +MagickSQ2;
1422
CalcKernelMetaData(kernel);
1423
ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1424
break;
1425
case 13:
1426
kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1427
if (kernel == (KernelInfo *) NULL)
1428
return(kernel);
1429
kernel->type = type;
1430
kernel->values[0] = +MagickSQ2;
1431
kernel->values[8] = -MagickSQ2;
1432
CalcKernelMetaData(kernel);
1433
ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1434
break;
1435
case 14:
1436
kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1437
if (kernel == (KernelInfo *) NULL)
1438
return(kernel);
1439
kernel->type = type;
1440
kernel->values[2] = -MagickSQ2;
1441
kernel->values[6] = +MagickSQ2;
1442
CalcKernelMetaData(kernel);
1443
ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1444
break;
1445
case 15:
1446
kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1447
if (kernel == (KernelInfo *) NULL)
1448
return(kernel);
1449
kernel->type = type;
1450
ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1451
break;
1452
case 16:
1453
kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1454
if (kernel == (KernelInfo *) NULL)
1455
return(kernel);
1456
kernel->type = type;
1457
ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1458
break;
1459
case 17:
1460
kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1461
if (kernel == (KernelInfo *) NULL)
1462
return(kernel);
1463
kernel->type = type;
1464
ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1465
break;
1466
case 18:
1467
kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1468
if (kernel == (KernelInfo *) NULL)
1469
return(kernel);
1470
kernel->type = type;
1471
ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1472
break;
1473
case 19:
1474
kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1475
if (kernel == (KernelInfo *) NULL)
1476
return(kernel);
1477
kernel->type = type;
1478
ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1479
break;
1480
}
1481
if ( fabs(args->sigma) > MagickEpsilon )
1482
/* Rotate by correctly supplied 'angle' */
1483
RotateKernelInfo(kernel, args->sigma);
1484
else if ( args->rho > 30.0 || args->rho < -30.0 )
1485
/* Rotate by out of bounds 'type' */
1486
RotateKernelInfo(kernel, args->rho);
1487
break;
1488
}
1489
1490
712
/* Boolean Kernels */
1491
case DiamondKernel:
713
case RectangleKernel:
714
case SquareKernel:
1492
715
{
1493
if (args->rho < 1.0)
1494
kernel->width = kernel->height = 3; /* default radius = 1 */
1495
else
1496
kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1497
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1498
1499
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1500
kernel->height*sizeof(double));
1501
if (kernel->values == (double *) NULL)
1502
return(DestroyKernelInfo(kernel));
1503
1504
/* set all kernel values within diamond area to scale given */
1505
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1506
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1507
if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1508
kernel->positive_range += kernel->values[i] = args->sigma;
1509
else
1510
kernel->values[i] = nan;
1511
kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1512
break;
1513
}
1514
case SquareKernel:
1515
case RectangleKernel:
1516
{ double
1517
scale;
716
double scale;
1518
717
if ( type == SquareKernel )
1519
718
{
1520
719
if (args->rho < 1.0)
1521
720
kernel->width = kernel->height = 3; /* default radius = 1 */
1522
721
else
1523
kernel->width = kernel->height = (size_t) (2*args->rho+1);
1524
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
722
kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
723
kernel->x = kernel->y = (long) (kernel->width-1)/2;
1525
724
scale = args->sigma;
1526
725
}
1527
726
else {
1528
727
/* NOTE: user defaults set in "AcquireKernelInfo()" */
1529
728
if ( args->rho < 1.0 || args->sigma < 1.0 )
1530
729
return(DestroyKernelInfo(kernel)); /* invalid args given */
1531
kernel->width = (size_t)args->rho;
1532
kernel->height = (size_t)args->sigma;
730
kernel->width = (unsigned long)args->rho;
731
kernel->height = (unsigned long)args->sigma;
1533
732
if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1534
733
args->psi < 0.0 || args->psi > (double)kernel->height )
1535
734
return(DestroyKernelInfo(kernel)); /* invalid args given */
1536
kernel->x = (ssize_t) args->xi;
1537
kernel->y = (ssize_t) args->psi;
735
kernel->x = (long) args->xi;
736
kernel->y = (long) args->psi;
1538
737
scale = 1.0;
1539
738
}
1540
739
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1542
741
if (kernel->values == (double *) NULL)
1543
742
return(DestroyKernelInfo(kernel));
1544
743
1545
/* set all kernel values to scale given */
1546
u=(ssize_t) (kernel->width*kernel->height);
744
/* set all kernel values to 1.0 */
745
u=(long) kernel->width*kernel->height;
1547
746
for ( i=0; i < u; i++)
1548
747
kernel->values[i] = scale;
1549
748
kernel->minimum = kernel->maximum = scale; /* a flat shape */
1550
749
kernel->positive_range = scale*u;
1551
750
break;
1552
751
}
752
case DiamondKernel:
753
{
754
if (args->rho < 1.0)
755
kernel->width = kernel->height = 3; /* default radius = 1 */
756
else
757
kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
758
kernel->x = kernel->y = (long) (kernel->width-1)/2;
759
760
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
761
kernel->height*sizeof(double));
762
if (kernel->values == (double *) NULL)
763
return(DestroyKernelInfo(kernel));
764
765
/* set all kernel values within diamond area to scale given */
766
for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
767
for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
768
if ((labs(u)+labs(v)) <= (long)kernel->x)
769
kernel->positive_range += kernel->values[i] = args->sigma;
770
else
771
kernel->values[i] = nan;
772
kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
773
break;
774
}
1553
775
case DiskKernel:
1554
776
{
1555
ssize_t
1556
limit = (ssize_t)(args->rho*args->rho);
777
long
778
limit;
1557
779
1558
if (args->rho < 0.4) /* default radius approx 3.5 */
780
limit = (long)(args->rho*args->rho);
781
if (args->rho < 0.1) /* default radius approx 3.5 */
1559
782
kernel->width = kernel->height = 7L, limit = 10L;
1560
783
else
1561
kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1562
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
784
kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
785
kernel->x = kernel->y = (long) (kernel->width-1)/2;
1563
786
1564
787
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1565
788
kernel->height*sizeof(double));
1566
789
if (kernel->values == (double *) NULL)
1567
790
return(DestroyKernelInfo(kernel));
1568
791
1569
/* set all kernel values within disk area to scale given */
1570
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1571
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
792
/* set all kernel values within disk area to 1.0 */
793
for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
794
for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1572
795
if ((u*u+v*v) <= limit)
1573
796
kernel->positive_range += kernel->values[i] = args->sigma;
1574
797
else
1581
804
if (args->rho < 1.0)
1582
805
kernel->width = kernel->height = 5; /* default radius 2 */
1583
806
else
1584
kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1585
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
807
kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
808
kernel->x = kernel->y = (long) (kernel->width-1)/2;
1586
809
1587
810
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1588
811
kernel->height*sizeof(double));
1589
812
if (kernel->values == (double *) NULL)
1590
813
return(DestroyKernelInfo(kernel));
1591
814
1592
/* set all kernel values along axises to given scale */
1593
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1594
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
815
/* set all kernel values along axises to 1.0 */
816
for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
817
for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1595
818
kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1596
819
kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1597
820
kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1598
821
break;
1599
822
}
1600
case CrossKernel:
1601
{
1602
if (args->rho < 1.0)
1603
kernel->width = kernel->height = 5; /* default radius 2 */
1604
else
1605
kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1606
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1607
1608
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1609
kernel->height*sizeof(double));
1610
if (kernel->values == (double *) NULL)
1611
return(DestroyKernelInfo(kernel));
1612
1613
/* set all kernel values along axises to given scale */
1614
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1615
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1616
kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1617
kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1618
kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1619
break;
1620
}
1621
/* HitAndMiss Kernels */
1622
case RingKernel:
1623
case PeaksKernel:
1624
{
1625
ssize_t
1626
limit1,
1627
limit2,
1628
scale;
1629
1630
if (args->rho < args->sigma)
1631
{
1632
kernel->width = ((size_t)args->sigma)*2+1;
1633
limit1 = (ssize_t)(args->rho*args->rho);
1634
limit2 = (ssize_t)(args->sigma*args->sigma);
1635
}
1636
else
1637
{
1638
kernel->width = ((size_t)args->rho)*2+1;
1639
limit1 = (ssize_t)(args->sigma*args->sigma);
1640
limit2 = (ssize_t)(args->rho*args->rho);
1641
}
1642
if ( limit2 <= 0 )
1643
kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1644
1645
kernel->height = kernel->width;
1646
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1647
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1648
kernel->height*sizeof(double));
1649
if (kernel->values == (double *) NULL)
1650
return(DestroyKernelInfo(kernel));
1651
1652
/* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1653
scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1654
for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1655
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1656
{ ssize_t radius=u*u+v*v;
1657
if (limit1 < radius && radius <= limit2)
1658
kernel->positive_range += kernel->values[i] = (double) scale;
1659
else
1660
kernel->values[i] = nan;
1661
}
1662
kernel->minimum = kernel->minimum = (double) scale;
1663
if ( type == PeaksKernel ) {
1664
/* set the central point in the middle */
1665
kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1666
kernel->positive_range = 1.0;
1667
kernel->maximum = 1.0;
1668
}
1669
break;
1670
}
1671
case EdgesKernel:
1672
{
1673
kernel=ParseKernelArray("3: 0,0,0 -,1,- 1,1,1");
1674
if (kernel == (KernelInfo *) NULL)
1675
return(kernel);
1676
kernel->type = type;
1677
ExpandMirrorKernelInfo(kernel); /* mirror expansion of other kernels */
1678
break;
1679
}
1680
case CornersKernel:
1681
{
1682
kernel=ParseKernelArray("3: 0,0,- 0,1,1 -,1,-");
1683
if (kernel == (KernelInfo *) NULL)
1684
return(kernel);
1685
kernel->type = type;
1686
ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1687
break;
1688
}
1689
case ThinDiagonalsKernel:
1690
{
1691
switch ( (int) args->rho ) {
1692
case 0:
1693
default:
1694
{ KernelInfo
1695
*new_kernel;
1696
kernel=ParseKernelArray("3: 0,0,0 0,1,1 1,1,-");
1697
if (kernel == (KernelInfo *) NULL)
1698
return(kernel);
1699
kernel->type = type;
1700
new_kernel=ParseKernelArray("3: 0,0,1 0,1,1 0,1,-");
1701
if (new_kernel == (KernelInfo *) NULL)
1702
return(DestroyKernelInfo(kernel));
1703
new_kernel->type = type;
1704
LastKernelInfo(kernel)->next = new_kernel;
1705
ExpandMirrorKernelInfo(kernel);
1706
break;
1707
}
1708
case 1:
1709
kernel=ParseKernelArray("3: 0,0,0 0,1,1 1,1,-");
1710
if (kernel == (KernelInfo *) NULL)
1711
return(kernel);
1712
kernel->type = type;
1713
RotateKernelInfo(kernel, args->sigma);
1714
break;
1715
case 2:
1716
kernel=ParseKernelArray("3: 0,0,1 0,1,1 0,1,-");
1717
if (kernel == (KernelInfo *) NULL)
1718
return(kernel);
1719
kernel->type = type;
1720
RotateKernelInfo(kernel, args->sigma);
1721
break;
1722
}
1723
break;
1724
}
1725
case LineEndsKernel:
1726
{ /* Kernels for finding the end of thin lines */
1727
switch ( (int) args->rho ) {
1728
case 0:
1729
default:
1730
/* set of kernels to find all end of lines */
1731
kernel=AcquireKernelInfo("LineEnds:1>;LineEnds:2>");
1732
if (kernel == (KernelInfo *) NULL)
1733
return(kernel);
1734
break;
1735
case 1:
1736
/* kernel for 4-connected line ends - no rotation */
1737
kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1738
if (kernel == (KernelInfo *) NULL)
1739
return(kernel);
1740
kernel->type = type;
1741
RotateKernelInfo(kernel, args->sigma);
1742
break;
1743
case 2:
1744
/* kernel to add for 8-connected lines - no rotation */
1745
kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1746
if (kernel == (KernelInfo *) NULL)
1747
return(kernel);
1748
kernel->type = type;
1749
RotateKernelInfo(kernel, args->sigma);
1750
break;
1751
case 3:
1752
/* kernel to add for orthogonal line ends - does not find corners */
1753
kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1754
if (kernel == (KernelInfo *) NULL)
1755
return(kernel);
1756
kernel->type = type;
1757
RotateKernelInfo(kernel, args->sigma);
1758
break;
1759
case 4:
1760
/* traditional line end - fails on last T end */
1761
kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1762
if (kernel == (KernelInfo *) NULL)
1763
return(kernel);
1764
kernel->type = type;
1765
RotateKernelInfo(kernel, args->sigma);
1766
break;
1767
}
1768
break;
1769
}
1770
case LineJunctionsKernel:
1771
{ /* kernels for finding the junctions of multiple lines */
1772
switch ( (int) args->rho ) {
1773
case 0:
1774
default:
1775
/* set of kernels to find all line junctions */
1776
kernel=AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>");
1777
if (kernel == (KernelInfo *) NULL)
1778
return(kernel);
1779
break;
1780
case 1:
1781
/* Y Junction */
1782
kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1783
if (kernel == (KernelInfo *) NULL)
1784
return(kernel);
1785
kernel->type = type;
1786
RotateKernelInfo(kernel, args->sigma);
1787
break;
1788
case 2:
1789
/* Diagonal T Junctions */
1790
kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1791
if (kernel == (KernelInfo *) NULL)
1792
return(kernel);
1793
kernel->type = type;
1794
RotateKernelInfo(kernel, args->sigma);
1795
break;
1796
case 3:
1797
/* Orthogonal T Junctions */
1798
kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1799
if (kernel == (KernelInfo *) NULL)
1800
return(kernel);
1801
kernel->type = type;
1802
RotateKernelInfo(kernel, args->sigma);
1803
break;
1804
case 4:
1805
/* Diagonal X Junctions */
1806
kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1807
if (kernel == (KernelInfo *) NULL)
1808
return(kernel);
1809
kernel->type = type;
1810
RotateKernelInfo(kernel, args->sigma);
1811
break;
1812
case 5:
1813
/* Orthogonal X Junctions - minimal diamond kernel */
1814
kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1815
if (kernel == (KernelInfo *) NULL)
1816
return(kernel);
1817
kernel->type = type;
1818
RotateKernelInfo(kernel, args->sigma);
1819
break;
1820
}
1821
break;
1822
}
1823
case RidgesKernel:
1824
{ /* Ridges - Ridge finding kernels */
1825
KernelInfo
1826
*new_kernel;
1827
switch ( (int) args->rho ) {
1828
case 1:
1829
default:
1830
kernel=ParseKernelArray("3x1:0,1,0");
1831
if (kernel == (KernelInfo *) NULL)
1832
return(kernel);
1833
kernel->type = type;
1834
ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1835
break;
1836
case 2:
1837
kernel=ParseKernelArray("4x1:0,1,1,0");
1838
if (kernel == (KernelInfo *) NULL)
1839
return(kernel);
1840
kernel->type = type;
1841
ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1842
1843
/* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1844
/* Unfortunatally we can not yet rotate a non-square kernel */
1845
/* But then we can't flip a non-symetrical kernel either */
1846
new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1847
if (new_kernel == (KernelInfo *) NULL)
1848
return(DestroyKernelInfo(kernel));
1849
new_kernel->type = type;
1850
LastKernelInfo(kernel)->next = new_kernel;
1851
new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1852
if (new_kernel == (KernelInfo *) NULL)
1853
return(DestroyKernelInfo(kernel));
1854
new_kernel->type = type;
1855
LastKernelInfo(kernel)->next = new_kernel;
1856
new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1857
if (new_kernel == (KernelInfo *) NULL)
1858
return(DestroyKernelInfo(kernel));
1859
new_kernel->type = type;
1860
LastKernelInfo(kernel)->next = new_kernel;
1861
new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1862
if (new_kernel == (KernelInfo *) NULL)
1863
return(DestroyKernelInfo(kernel));
1864
new_kernel->type = type;
1865
LastKernelInfo(kernel)->next = new_kernel;
1866
new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1867
if (new_kernel == (KernelInfo *) NULL)
1868
return(DestroyKernelInfo(kernel));
1869
new_kernel->type = type;
1870
LastKernelInfo(kernel)->next = new_kernel;
1871
new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1872
if (new_kernel == (KernelInfo *) NULL)
1873
return(DestroyKernelInfo(kernel));
1874
new_kernel->type = type;
1875
LastKernelInfo(kernel)->next = new_kernel;
1876
new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1877
if (new_kernel == (KernelInfo *) NULL)
1878
return(DestroyKernelInfo(kernel));
1879
new_kernel->type = type;
1880
LastKernelInfo(kernel)->next = new_kernel;
1881
new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1882
if (new_kernel == (KernelInfo *) NULL)
1883
return(DestroyKernelInfo(kernel));
1884
new_kernel->type = type;
1885
LastKernelInfo(kernel)->next = new_kernel;
1886
break;
1887
}
1888
break;
1889
}
1890
case ConvexHullKernel:
1891
{
1892
KernelInfo
1893
*new_kernel;
1894
/* first set of 8 kernels */
1895
kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1896
if (kernel == (KernelInfo *) NULL)
1897
return(kernel);
1898
kernel->type = type;
1899
ExpandRotateKernelInfo(kernel, 90.0);
1900
/* append the mirror versions too - no flip function yet */
1901
new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1902
if (new_kernel == (KernelInfo *) NULL)
1903
return(DestroyKernelInfo(kernel));
1904
new_kernel->type = type;
1905
ExpandRotateKernelInfo(new_kernel, 90.0);
1906
LastKernelInfo(kernel)->next = new_kernel;
1907
break;
1908
}
1909
case SkeletonKernel:
1910
{
1911
KernelInfo
1912
*new_kernel;
1913
switch ( (int) args->rho ) {
1914
case 1:
1915
default:
1916
/* Traditional Skeleton...
1917
** A cyclically rotated single kernel
1918
*/
1919
kernel=ParseKernelArray("3: 0,0,0 -,1,- 1,1,1");
1920
if (kernel == (KernelInfo *) NULL)
1921
return(kernel);
1922
kernel->type = type;
1923
ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1924
break;
1925
case 2:
1926
/* HIPR Variation of the cyclic skeleton
1927
** Corners of the traditional method made more forgiving,
1928
** but the retain the same cyclic order.
1929
*/
1930
kernel=ParseKernelArray("3: 0,0,0 -,1,- 1,1,1");
1931
if (kernel == (KernelInfo *) NULL)
1932
return(kernel);
1933
kernel->type = type;
1934
new_kernel=ParseKernelArray("3: -,0,0 1,1,0 -,1,-");
1935
if (new_kernel == (KernelInfo *) NULL)
1936
return(new_kernel);
1937
new_kernel->type = type;
1938
LastKernelInfo(kernel)->next = new_kernel;
1939
ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1940
break;
1941
case 3:
1942
/* Jittered Skeleton: do top, then bottom, then each sides */
1943
/* Do top edge */
1944
kernel=ParseKernelArray("3: 0,0,0 -,1,- 1,1,1");
1945
if (kernel == (KernelInfo *) NULL)
1946
return(kernel);
1947
kernel->type = type;
1948
new_kernel=ParseKernelArray("3: 0,0,- 0,1,1 -,1,-");
1949
if (new_kernel == (KernelInfo *) NULL)
1950
return(new_kernel);
1951
new_kernel->type = type;
1952
LastKernelInfo(kernel)->next = new_kernel;
1953
new_kernel=ParseKernelArray("3: -,0,0 1,1,0 -,1,-");
1954
if (new_kernel == (KernelInfo *) NULL)
1955
return(new_kernel);
1956
new_kernel->type = type;
1957
LastKernelInfo(kernel)->next = new_kernel;
1958
/* Do Bottom edge */
1959
new_kernel=ParseKernelArray("3: 1,1,1 -,1,- 0,0,0");
1960
if (new_kernel == (KernelInfo *) NULL)
1961
return(new_kernel);
1962
new_kernel->type = type;
1963
LastKernelInfo(kernel)->next = new_kernel;
1964
new_kernel=ParseKernelArray("3: -,1,- 1,1,0 -,0,0");
1965
if (new_kernel == (KernelInfo *) NULL)
1966
return(new_kernel);
1967
new_kernel->type = type;
1968
LastKernelInfo(kernel)->next = new_kernel;
1969
new_kernel=ParseKernelArray("3: -,1,- 0,1,1 0,0,-");
1970
if (new_kernel == (KernelInfo *) NULL)
1971
return(new_kernel);
1972
new_kernel->type = type;
1973
LastKernelInfo(kernel)->next = new_kernel;
1974
/* Last the two sides */
1975
new_kernel=ParseKernelArray("3: 0,-,1 0,1,1 0,-,1");
1976
if (new_kernel == (KernelInfo *) NULL)
1977
return(new_kernel);
1978
new_kernel->type = type;
1979
LastKernelInfo(kernel)->next = new_kernel;
1980
new_kernel=ParseKernelArray("3: 1,-,0 1,1,0 1,-,0");
1981
if (new_kernel == (KernelInfo *) NULL)
1982
return(new_kernel);
1983
new_kernel->type = type;
1984
LastKernelInfo(kernel)->next = new_kernel;
1985
break;
1986
case 4:
1987
/* Just a simple 'Edge' kernel, but with a extra two kernels
1988
** to finish off diagonal lines, top then bottom then sides.
1989
** Works well for test case but fails for general case.
1990
*/
1991
kernel=ParseKernelArray("3: 0,0,0 -,1,- 1,1,1");
1992
if (kernel == (KernelInfo *) NULL)
1993
return(kernel);
1994
kernel->type = type;
1995
new_kernel=ParseKernelArray("3: 0,0,0 0,1,1 1,1,-");
1996
if (new_kernel == (KernelInfo *) NULL)
1997
return(DestroyKernelInfo(kernel));
1998
new_kernel->type = type;
1999
LastKernelInfo(kernel)->next = new_kernel;
2000
new_kernel=ParseKernelArray("3: 0,0,0 1,1,0 -,1,1");
2001
if (new_kernel == (KernelInfo *) NULL)
2002
return(DestroyKernelInfo(kernel));
2003
new_kernel->type = type;
2004
LastKernelInfo(kernel)->next = new_kernel;
2005
ExpandMirrorKernelInfo(kernel);
2006
/* Append a set of corner kernels */
2007
new_kernel=ParseKernelArray("3: 0,0,- 0,1,1 -,1,-");
2008
if (new_kernel == (KernelInfo *) NULL)
2009
return(DestroyKernelInfo(kernel));
2010
new_kernel->type = type;
2011
ExpandRotateKernelInfo(new_kernel, 90.0);
2012
LastKernelInfo(kernel)->next = new_kernel;
2013
break;
2014
}
2015
break;
2016
}
2017
823
/* Distance Measuring Kernels */
2018
824
case ChebyshevKernel:
2019
825
{
826
double
827
scale;
828
2020
829
if (args->rho < 1.0)
2021
830
kernel->width = kernel->height = 3; /* default radius = 1 */
2022
831
else
2023
kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2024
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
832
kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
833
kernel->x = kernel->y = (long) (kernel->width-1)/2;
2025
834
2026
835
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2027
836
kernel->height*sizeof(double));
2028
837
if (kernel->values == (double *) NULL)
2029
838
return(DestroyKernelInfo(kernel));
2030
839
2031
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2032
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
840
scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
841
for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
842
for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
2033
843
kernel->positive_range += ( kernel->values[i] =
2034
args->sigma*((labs((long) u)>labs((long) v)) ? labs((long) u) : labs((long) v)) );
844
scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
2035
845
kernel->maximum = kernel->values[0];
2036
846
break;
2037
847
}
2038
case ManhattanKernel:
848
case ManhattenKernel:
2039
849
{
850
double
851
scale;
852
2040
853
if (args->rho < 1.0)
2041
854
kernel->width = kernel->height = 3; /* default radius = 1 */
2042
855
else
2043
kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2044
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
856
kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
857
kernel->x = kernel->y = (long) (kernel->width-1)/2;
2045
858
2046
859
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2047
860
kernel->height*sizeof(double));
2048
861
if (kernel->values == (double *) NULL)
2049
862
return(DestroyKernelInfo(kernel));
2050
863
2051
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2052
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
864
scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
865
for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
866
for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
2053
867
kernel->positive_range += ( kernel->values[i] =
2054
args->sigma*(labs((long) u)+labs((long) v)) );
868
scale*(labs(u)+labs(v)) );
2055
869
kernel->maximum = kernel->values[0];
2056
870
break;
2057
871
}
2058
872
case EuclideanKernel:
2059
873
{
874
double
875
scale;
876
2060
877
if (args->rho < 1.0)
2061
878
kernel->width = kernel->height = 3; /* default radius = 1 */
2062
879
else
2063
kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2064
kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
880
kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
881
kernel->x = kernel->y = (long) (kernel->width-1)/2;
2065
882
2066
883
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2067
884
kernel->height*sizeof(double));
2068
885
if (kernel->values == (double *) NULL)
2069
886
return(DestroyKernelInfo(kernel));
2070
887
2071
for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2072
for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
888
scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
889
for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
890
for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
2073
891
kernel->positive_range += ( kernel->values[i] =
2074
args->sigma*sqrt((double)(u*u+v*v)) );
892
scale*sqrt((double)(u*u+v*v)) );
2075
893
kernel->maximum = kernel->values[0];
2076
894
break;
2077
895
}
2078
case UnityKernel:
896
/* Undefined Kernels */
897
case LaplacianKernel:
898
case LOGKernel:
899
case DOGKernel:
900
perror("Kernel Type has not been defined yet");
901
/* FALL THRU */
2079
902
default:
2080
{
2081
/* Unity or No-Op Kernel - Basically just a single pixel on its own */
2082
kernel=ParseKernelArray("1:1");
2083
if (kernel == (KernelInfo *) NULL)
2084
return(kernel);
2085
kernel->type = ( type == UnityKernel ) ? UnityKernel : UndefinedKernel;
2086
break;
2087
}
903
/* Generate a No-Op minimal kernel - 1x1 pixel */
904
kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
905
if (kernel->values == (double *) NULL)
906
return(DestroyKernelInfo(kernel));
907
kernel->width = kernel->height = 1;
908
kernel->x = kernel->x = 0;
909
kernel->type = UndefinedKernel;
910
kernel->maximum =
911
kernel->positive_range =
912
kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
2088
913
break;
2089
914
}
2090
915
2102
927
% %
2103
928
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2104
929
%
2105
% CloneKernelInfo() creates a new clone of the given Kernel List so that its
2106
% can be modified without effecting the original. The cloned kernel should
2107
% be destroyed using DestoryKernelInfo() when no ssize_ter needed.
930
% CloneKernelInfo() creates a new clone of the given Kernel so that its can
931
% be modified without effecting the original. The cloned kernel should be
932
% destroyed using DestoryKernelInfo() when no longer needed.
2108
933
%
2109
% The format of the CloneKernelInfo method is:
934
% The format of the DestroyKernelInfo method is:
2110
935
%
2111
936
% KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2112
937
%
2117
942
*/
2118
943
MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2119
944
{
2120
register ssize_t
945
register long
2121
946
i;
2122
947
2123
KernelInfo
2124
*new_kernel;
948
KernelInfo *
949
new;
2125
950
2126
951
assert(kernel != (KernelInfo *) NULL);
2127
new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2128
if (new_kernel == (KernelInfo *) NULL)
2129
return(new_kernel);
2130
*new_kernel=(*kernel); /* copy values in structure */
2131
2132
/* replace the values with a copy of the values */
2133
new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2134
kernel->height*sizeof(double));
2135
if (new_kernel->values == (double *) NULL)
2136
return(DestroyKernelInfo(new_kernel));
2137
for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2138
new_kernel->values[i]=kernel->values[i];
2139
2140
/* Also clone the next kernel in the kernel list */
2141
if ( kernel->next != (KernelInfo *) NULL ) {
2142
new_kernel->next = CloneKernelInfo(kernel->next);
2143
if ( new_kernel->next == (KernelInfo *) NULL )
2144
return(DestroyKernelInfo(new_kernel));
2145
}
2146
2147
return(new_kernel);
952
953
new=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
954
if (new == (KernelInfo *) NULL)
955
return(new);
956
*new = *kernel; /* copy values in structure */
957
958
new->values=(double *) AcquireQuantumMemory(kernel->width,
959
kernel->height*sizeof(double));
960
if (new->values == (double *) NULL)
961
return(DestroyKernelInfo(new));
962
963
for (i=0; i < (long) (kernel->width*kernel->height); i++)
964
new->values[i] = kernel->values[i];
965
966
return(new);
2148
967
}
2149
968
2150
969
/*
2170
989
% o kernel: the Morphology/Convolution kernel to be destroyed
2171
990
%
2172
991
*/
992
2173
993
MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2174
994
{
2175
995
assert(kernel != (KernelInfo *) NULL);
2176
996
2177
if ( kernel->next != (KernelInfo *) NULL )
2178
kernel->next = DestroyKernelInfo(kernel->next);
2179
2180
kernel->values = (double *)RelinquishMagickMemory(kernel->values);
2181
kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
997
kernel->values=(double *) AcquireQuantumMemory(kernel->width,
998
kernel->height*sizeof(double));
999
kernel->values=(double *)RelinquishMagickMemory(kernel->values);
1000
kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2182
1001
return(kernel);
2183
1002
}
2184
1003
2187
1006
% %
2188
1007
% %
2189
1008
% %
2190
% E x p a n d M i r r o r K e r n e l I n f o %
2191
% %
2192
% %
2193
% %
2194
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2195
%
2196
% ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2197
% sequence of 90-degree rotated kernels but providing a reflected 180
2198
% rotatation, before the -/+ 90-degree rotations.
2199
%
2200
% This special rotation order produces a better, more symetrical thinning of
2201
% objects.
2202
%
2203
% The format of the ExpandMirrorKernelInfo method is:
2204
%
2205
% void ExpandMirrorKernelInfo(KernelInfo *kernel)
2206
%
2207
% A description of each parameter follows:
2208
%
2209
% o kernel: the Morphology/Convolution kernel
2210
%
2211
% This function is only internel to this module, as it is not finalized,
2212
% especially with regard to non-orthogonal angles, and rotation of larger
2213
% 2D kernels.
2214
*/
2215
2216
#if 0
2217
static void FlopKernelInfo(KernelInfo *kernel)
2218
{ /* Do a Flop by reversing each row. */
2219
size_t
2220
y;
2221
register ssize_t
2222
x,r;
2223
register double
2224
*k,t;
2225
2226
for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2227
for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2228
t=k[x], k[x]=k[r], k[r]=t;
2229
2230
kernel->x = kernel->width - kernel->x - 1;
2231
angle = fmod(angle+180.0, 360.0);
2232
}
2233
#endif
2234
2235
static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2236
{
2237
KernelInfo
2238
*clone,
2239
*last;
2240
2241
last = kernel;
2242
2243
clone = CloneKernelInfo(last);
2244
RotateKernelInfo(clone, 180); /* flip */
2245
LastKernelInfo(last)->next = clone;
2246
last = clone;
2247
2248
clone = CloneKernelInfo(last);
2249
RotateKernelInfo(clone, 90); /* transpose */
2250
LastKernelInfo(last)->next = clone;
2251
last = clone;
2252
2253
clone = CloneKernelInfo(last);
2254
RotateKernelInfo(clone, 180); /* flop */
2255
LastKernelInfo(last)->next = clone;
2256
2257
return;
2258
}
2259
2260
/*
2261
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2262
% %
2263
% %
2264
% %
2265
% E x p a n d R o t a t e K e r n e l I n f o %
2266
% %
2267
% %
2268
% %
2269
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2270
%
2271
% ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2272
% incrementally by the angle given, until the first kernel repeats.
2273
%
2274
% WARNING: 45 degree rotations only works for 3x3 kernels.
2275
% While 90 degree roatations only works for linear and square kernels
2276
%
2277
% The format of the ExpandRotateKernelInfo method is:
2278
%
2279
% void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2280
%
2281
% A description of each parameter follows:
2282
%
2283
% o kernel: the Morphology/Convolution kernel
2284
%
2285
% o angle: angle to rotate in degrees
2286
%
2287
% This function is only internel to this module, as it is not finalized,
2288
% especially with regard to non-orthogonal angles, and rotation of larger
2289
% 2D kernels.
2290
*/
2291
2292
/* Internal Routine - Return true if two kernels are the same */
2293
static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2294
const KernelInfo *kernel2)
2295
{
2296
register size_t
2297
i;
2298
2299
/* check size and origin location */
2300
if ( kernel1->width != kernel2->width
2301
|| kernel1->height != kernel2->height
2302
|| kernel1->x != kernel2->x
2303
|| kernel1->y != kernel2->y )
2304
return MagickFalse;
2305
2306
/* check actual kernel values */
2307
for (i=0; i < (kernel1->width*kernel1->height); i++) {
2308
/* Test for Nan equivelence */
2309
if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2310
return MagickFalse;
2311
if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2312
return MagickFalse;
2313
/* Test actual values are equivelent */
2314
if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
2315
return MagickFalse;
2316
}
2317
2318
return MagickTrue;
2319
}
2320
2321
static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2322
{
2323
KernelInfo
2324
*clone,
2325
*last;
2326
2327
last = kernel;
2328
while(1) {
2329
clone = CloneKernelInfo(last);
2330
RotateKernelInfo(clone, angle);
2331
if ( SameKernelInfo(kernel, clone) == MagickTrue )
2332
break;
2333
LastKernelInfo(last)->next = clone;
2334
last = clone;
2335
}
2336
clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2337
return;
2338
}
2339
2340
/*
2341
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2342
% %
2343
% %
2344
% %
2345
+ C a l c M e t a K e r n a l I n f o %
2346
% %
2347
% %
2348
% %
2349
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2350
%
2351
% CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2352
% using the kernel values. This should only ne used if it is not posible to
2353
% calculate that meta-data in some easier way.
2354
%
2355
% It is important that the meta-data is correct before ScaleKernelInfo() is
2356
% used to perform kernel normalization.
2357
%
2358
% The format of the CalcKernelMetaData method is:
2359
%
2360
% void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2361
%
2362
% A description of each parameter follows:
2363
%
2364
% o kernel: the Morphology/Convolution kernel to modify
2365
%
2366
% WARNING: Minimum and Maximum values are assumed to include zero, even if
2367
% zero is not part of the kernel (as in Gaussian Derived kernels). This
2368
% however is not true for flat-shaped morphological kernels.
2369
%
2370
% WARNING: Only the specific kernel pointed to is modified, not a list of
2371
% multiple kernels.
2372
%
2373
% This is an internal function and not expected to be useful outside this
2374
% module. This could change however.
2375
*/
2376
static void CalcKernelMetaData(KernelInfo *kernel)
2377
{
2378
register size_t
2379
i;
2380
2381
kernel->minimum = kernel->maximum = 0.0;
2382
kernel->negative_range = kernel->positive_range = 0.0;
2383
for (i=0; i < (kernel->width*kernel->height); i++)
2384
{
2385
if ( fabs(kernel->values[i]) < MagickEpsilon )
2386
kernel->values[i] = 0.0;
2387
( kernel->values[i] < 0)
2388
? ( kernel->negative_range += kernel->values[i] )
2389
: ( kernel->positive_range += kernel->values[i] );
2390
Minimize(kernel->minimum, kernel->values[i]);
2391
Maximize(kernel->maximum, kernel->values[i]);
2392
}
2393
2394
return;
2395
}
2396
2397
/*
2398
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2399
% %
2400
% %
2401
% %
2402
% M o r p h o l o g y A p p l y %
2403
% %
2404
% %
2405
% %
2406
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2407
%
2408
% MorphologyApply() applies a morphological method, multiple times using
2409
% a list of multiple kernels.
2410
%
2411
% It is basically equivelent to as MorphologyImageChannel() (see below) but
2412
% without any user controls. This allows internel programs to use this
2413
% function, to actually perform a specific task without posible interference
2414
% by any API user supplied settings.
2415
%
2416
% It is MorphologyImageChannel() task to extract any such user controls, and
2417
% pass them to this function for processing.
2418
%
2419
% More specifically kernels are not normalized/scaled/blended by the
2420
% 'convolve:scale' Image Artifact (setting), nor is the convolve bias
2421
% (-bias setting or image->bias) loooked at, but must be supplied from the
2422
% function arguments.
2423
%
2424
% The format of the MorphologyApply method is:
2425
%
2426
% Image *MorphologyApply(const Image *image,MorphologyMethod method,
2427
% const ssize_t iterations,const KernelInfo *kernel,
2428
% const CompositeMethod compose, const double bias,
2429
% ExceptionInfo *exception)
2430
%
2431
% A description of each parameter follows:
2432
%
2433
% o image: the source image
1009
% M o r p h o l o g y I m a g e C h a n n e l %
1010
% %
1011
% %
1012
% %
1013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1014
%
1015
% MorphologyImageChannel() applies a user supplied kernel to the image
1016
% according to the given mophology method.
1017
%
1018
% The given kernel is assumed to have been pre-scaled appropriatally, usally
1019
% by the kernel generator.
1020
%
1021
% The format of the MorphologyImage method is:
1022
%
1023
% Image *MorphologyImage(const Image *image,MorphologyMethod method,
1024
% const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
1025
% Image *MorphologyImageChannel(const Image *image, const ChannelType
1026
% channel,MorphologyMethod method,const long iterations,
1027
% KernelInfo *kernel,ExceptionInfo *exception)
1028
%
1029
% A description of each parameter follows:
1030
%
1031
% o image: the image.
2434
1032
%
2435
1033
% o method: the morphology method to be applied.
2436
1034
%
2442
1040
% o channel: the channel type.
2443
1041
%
2444
1042
% o kernel: An array of double representing the morphology kernel.
2445
%
2446
% o compose: How to handle or merge multi-kernel results.
2447
% If 'UndefinedCompositeOp' use default for the Morphology method.
2448
% If 'NoCompositeOp' force image to be re-iterated by each kernel.
2449
% Otherwise merge the results using the compose method given.
2450
%
2451
% o bias: Convolution Output Bias.
1043
% Warning: kernel may be normalized for the Convolve method.
2452
1044
%
2453
1045
% o exception: return any errors or warnings in this structure.
2454
1046
%
2455
*/
2456
2457
2458
/* Apply a Morphology Primative to an image using the given kernel.
2459
** Two pre-created images must be provided, no image is created.
2460
** It returns the number of pixels that changed betwene the images
2461
** for convergence determination.
2462
*/
2463
static size_t MorphologyPrimitive(const Image *image, Image
1047
%
1048
% TODO: bias and auto-scale handling of the kernel for convolution
1049
% The given kernel is assumed to have been pre-scaled appropriatally, usally
1050
% by the kernel generator.
1051
%
1052
*/
1053
1054
1055
/* Internal function
1056
* Apply the Low-Level Morphology Method using the given Kernel
1057
* Returning the number of pixels that changed.
1058
* Two pre-created images must be provided, no image is created.
1059
*/
1060
static unsigned long MorphologyApply(const Image *image, Image
2464
1061
*result_image, const MorphologyMethod method, const ChannelType channel,
2465
const KernelInfo *kernel,const double bias,ExceptionInfo *exception)
1062
const KernelInfo *kernel, ExceptionInfo *exception)
2466
1063
{
2467
1064
#define MorphologyTag "Morphology/Image"
2468
1065
1066
long
1067
progress,
1068
y, offx, offy,
1069
changed;
1070
1071
MagickBooleanType
1072
status;
1073
1074
MagickPixelPacket
1075
bias;
1076
2469
1077
CacheView
2470
1078
*p_view,
2471
1079
*q_view;
2472
1080
2473
ssize_t
2474
y, offx, offy,
2475
changed;
2476
2477
MagickBooleanType
2478
status;
2479
2480
MagickOffsetType
2481
progress;
2482
2483
assert(image != (Image *) NULL);
2484
assert(image->signature == MagickSignature);
2485
assert(result_image != (Image *) NULL);
2486
assert(result_image->signature == MagickSignature);
2487
assert(kernel != (KernelInfo *) NULL);
2488
assert(kernel->signature == MagickSignature);
2489
assert(exception != (ExceptionInfo *) NULL);
2490
assert(exception->signature == MagickSignature);
2491
1081
/* Only the most basic morphology is actually performed by this routine */
1082
1083
/*
1084
Apply Basic Morphology to Image.
1085
*/
2492
1086
status=MagickTrue;
2493
1087
changed=0;
2494
1088
progress=0;
2495
1089
1090
GetMagickPixelPacket(image,&bias);
1091
SetMagickPixelPacketBias(image,&bias);
1092
/* Future: handle auto-bias from user, based on kernel input */
1093
2496
1094
p_view=AcquireCacheView(image);
2497
1095
q_view=AcquireCacheView(result_image);
2498
1096
2499
1097
/* Some methods (including convolve) needs use a reflected kernel.
2500
* Adjust 'origin' offsets to loop though kernel as a reflection.
1098
* Adjust 'origin' offsets for this reflected kernel.
2501
1099
*/
2502
1100
offx = kernel->x;
2503
1101
offy = kernel->y;
2504
1102
switch(method) {
1103
case ErodeMorphology:
1104
case ErodeIntensityMorphology:
1105
/* kernel is user as is, without reflection */
1106
break;
2505
1107
case ConvolveMorphology:
2506
1108
case DilateMorphology:
2507
1109
case DilateIntensityMorphology:
2508
1110
case DistanceMorphology:
2509
/* kernel needs to used with reflection about origin */
2510
offx = (ssize_t) kernel->width-offx-1;
2511
offy = (ssize_t) kernel->height-offy-1;
2512
break;
2513
case ErodeMorphology:
2514
case ErodeIntensityMorphology:
2515
case HitAndMissMorphology:
2516
case ThinningMorphology:
2517
case ThickenMorphology:
2518
/* kernel is used as is, without reflection */
1111
/* kernel needs to used with reflection */
1112
offx = (long) kernel->width-offx-1;
1113
offy = (long) kernel->height-offy-1;
2519
1114
break;
2520
1115
default:
2521
assert("Not a Primitive Morphology Method" != (char *) NULL);
1116
perror("Not a low level Morpholgy Method");
2522
1117
break;
2523
1118
}
2524
1119
2525
2526
if ( method == ConvolveMorphology && kernel->width == 1 )
2527
{ /* Special handling (for speed) of vertical (blur) kernels.
2528
** This performs its handling in columns rather than in rows.
2529
** This is only done fo convolve as it is the only method that
2530
** generates very large 1-D vertical kernels (such as a 'BlurKernel')
2531
**
2532
** Timing tests (on single CPU laptop)
2533
** Using a vertical 1-d Blue with normal row-by-row (below)
2534
** time convert logo: -morphology Convolve Blur:0x10+90 null:
2535
** 0.807u
2536
** Using this column method
2537
** time convert logo: -morphology Convolve Blur:0x10+90 null:
2538
** 0.620u
2539
**
2540
** Anthony Thyssen, 14 June 2010
2541
*/
2542
register ssize_t
2543
x;
2544
2545
#if defined(MAGICKCORE_OPENMP_SUPPORT)
2546
#pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2547
#endif
2548
for (x=0; x < (ssize_t) image->columns; x++)
2549
{
2550
register const PixelPacket
2551
*restrict p;
2552
2553
register const IndexPacket
2554
*restrict p_indexes;
2555
2556
register PixelPacket
2557
*restrict q;
2558
2559
register IndexPacket
2560
*restrict q_indexes;
2561
2562
register ssize_t
2563
y;
2564
2565
size_t
2566
r;
2567
2568
if (status == MagickFalse)
2569
continue;
2570
p=GetCacheViewVirtualPixels(p_view, x, -offy,1,
2571
image->rows+kernel->height, exception);
2572
q=GetCacheViewAuthenticPixels(q_view,x,0,1,result_image->rows,exception);
2573
if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2574
{
2575
status=MagickFalse;
2576
continue;
2577
}
2578
p_indexes=GetCacheViewVirtualIndexQueue(p_view);
2579
q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
2580
r = offy; /* offset to the origin pixel in 'p' */
2581
2582
for (y=0; y < (ssize_t) image->rows; y++)
2583
{
2584
register ssize_t
2585
v;
2586
2587
register const double
2588
*restrict k;
2589
2590
register const PixelPacket
2591
*restrict k_pixels;
2592
2593
register const IndexPacket
2594
*restrict k_indexes;
2595
2596
MagickPixelPacket
2597
result;
2598
2599
/* Copy input image to the output image for unused channels
2600
* This removes need for 'cloning' a new image every iteration
2601
*/
2602
*q = p[r];
2603
if (image->colorspace == CMYKColorspace)
2604
q_indexes[y] = p_indexes[r];
2605
2606
/* Set the bias of the weighted average output */
2607
result.red =
2608
result.green =
2609
result.blue =
2610
result.opacity =
2611
result.index = bias;
2612
2613
2614
/* Weighted Average of pixels using reflected kernel
2615
**
2616
** NOTE for correct working of this operation for asymetrical
2617
** kernels, the kernel needs to be applied in its reflected form.
2618
** That is its values needs to be reversed.
2619
*/
2620
k = &kernel->values[ kernel->height-1 ];
2621
k_pixels = p;
2622
k_indexes = p_indexes;
2623
if ( ((channel & SyncChannels) == 0 ) ||
2624
(image->matte == MagickFalse) )
2625
{ /* No 'Sync' involved.
2626
** Convolution is simple greyscale channel operation
2627
*/
2628
for (v=0; v < (ssize_t) kernel->height; v++) {
2629
if ( IsNan(*k) ) continue;
2630
result.red += (*k)*k_pixels->red;
2631
result.green += (*k)*k_pixels->green;
2632
result.blue += (*k)*k_pixels->blue;
2633
result.opacity += (*k)*k_pixels->opacity;
2634
if ( image->colorspace == CMYKColorspace)
2635
result.index += (*k)*(*k_indexes);
2636
k--;
2637
k_pixels++;
2638
k_indexes++;
2639
}
2640
if ((channel & RedChannel) != 0)
2641
q->red = ClampToQuantum(result.red);
2642
if ((channel & GreenChannel) != 0)
2643
q->green = ClampToQuantum(result.green);
2644
if ((channel & BlueChannel) != 0)
2645
q->blue = ClampToQuantum(result.blue);
2646
if ((channel & OpacityChannel) != 0
2647
&& image->matte == MagickTrue )
2648
q->opacity = ClampToQuantum(result.opacity);
2649
if ((channel & IndexChannel) != 0
2650
&& image->colorspace == CMYKColorspace)
2651
q_indexes[x] = ClampToQuantum(result.index);
2652
}
2653
else
2654
{ /* Channel 'Sync' Flag, and Alpha Channel enabled.
2655
** Weight the color channels with Alpha Channel so that
2656
** transparent pixels are not part of the results.
2657
*/
2658
MagickRealType
2659
alpha, /* alpha weighting of colors : kernel*alpha */
2660
gamma; /* divisor, sum of color weighting values */
2661
2662
gamma=0.0;
2663
for (v=0; v < (ssize_t) kernel->height; v++) {
2664
if ( IsNan(*k) ) continue;
2665
alpha=(*k)*(QuantumScale*(QuantumRange-k_pixels->opacity));
2666
gamma += alpha;
2667
result.red += alpha*k_pixels->red;
2668
result.green += alpha*k_pixels->green;
2669
result.blue += alpha*k_pixels->blue;
2670
result.opacity += (*k)*k_pixels->opacity;
2671
if ( image->colorspace == CMYKColorspace)
2672
result.index += alpha*(*k_indexes);
2673
k--;
2674
k_pixels++;
2675
k_indexes++;
2676
}
2677
/* Sync'ed channels, all channels are modified */
2678
gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2679
q->red = ClampToQuantum(gamma*result.red);
2680
q->green = ClampToQuantum(gamma*result.green);
2681
q->blue = ClampToQuantum(gamma*result.blue);
2682
q->opacity = ClampToQuantum(result.opacity);
2683
if (image->colorspace == CMYKColorspace)
2684
q_indexes[x] = ClampToQuantum(gamma*result.index);
2685
}
2686
2687
/* Count up changed pixels */
2688
if ( ( p[r].red != q->red )
2689
|| ( p[r].green != q->green )
2690
|| ( p[r].blue != q->blue )
2691
|| ( p[r].opacity != q->opacity )
2692
|| ( image->colorspace == CMYKColorspace &&
2693
p_indexes[r] != q_indexes[x] ) )
2694
changed++; /* The pixel was changed in some way! */
2695
p++;
2696
q++;
2697
} /* y */
2698
if ( SyncCacheViewAuthenticPixels(q_view,exception) == MagickFalse)
2699
status=MagickFalse;
2700
if (image->progress_monitor != (MagickProgressMonitor) NULL)
2701
{
2702
MagickBooleanType
2703
proceed;
2704
2705
#if defined(MAGICKCORE_OPENMP_SUPPORT)
2706
#pragma omp critical (MagickCore_MorphologyImage)
2707
#endif
2708
proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2709
if (proceed == MagickFalse)
2710
status=MagickFalse;
2711
}
2712
} /* x */
2713
result_image->type=image->type;
2714
q_view=DestroyCacheView(q_view);
2715
p_view=DestroyCacheView(p_view);
2716
return(status ? (size_t) changed : 0);
2717
}
2718
2719
/*
2720
** Normal handling of horizontal or rectangular kernels (row by row)
2721
*/
2722
1120
#if defined(MAGICKCORE_OPENMP_SUPPORT)
2723
1121
#pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2724
1122
#endif
2725
for (y=0; y < (ssize_t) image->rows; y++)
1123
for (y=0; y < (long) image->rows; y++)
2726
1124
{
1125
MagickBooleanType
1126
sync;
1127
2727
1128
register const PixelPacket
2728
1129
*restrict p;
2729
1130
2736
1137
register IndexPacket
2737
1138
*restrict q_indexes;
2738
1139
2739
register ssize_t
1140
register long
2740
1141
x;
2741
1142
2742
size_t
1143
unsigned long
2743
1144
r;
2744
1145
2745
1146
if (status == MagickFalse)
2755
1156
}
2756
1157
p_indexes=GetCacheViewVirtualIndexQueue(p_view);
2757
1158
q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
2758
r = (image->columns+kernel->width)*offy+offx; /* offset to origin in 'p' */
1159
r = (image->columns+kernel->width)*offy+offx; /* constant */
2759
1160
2760
for (x=0; x < (ssize_t) image->columns; x++)
1161
for (x=0; x < (long) image->columns; x++)
2761
1162
{
2762
ssize_t
1163
long
2763
1164
v;
2764
1165
2765
register ssize_t
1166
register long
2766
1167
u;
2767
1168
2768
1169
register const double
2775
1176
*restrict k_indexes;
2776
1177
2777
1178
MagickPixelPacket
2778
result,
2779
min,
2780
max;
1179
result;
2781
1180
2782
/* Copy input image to the output image for unused channels
1181
/* Copy input to ouput image for unused channels
2783
1182
* This removes need for 'cloning' a new image every iteration
2784
1183
*/
2785
1184
*q = p[r];
2786
1185
if (image->colorspace == CMYKColorspace)
2787
1186
q_indexes[x] = p_indexes[r];
2788
1187
2789
/* Defaults */
2790
min.red =
2791
min.green =
2792
min.blue =
2793
min.opacity =
2794
min.index = (MagickRealType) QuantumRange;
2795
max.red =
2796
max.green =
2797
max.blue =
2798
max.opacity =
2799
max.index = (MagickRealType) 0;
2800
/* default result is the original pixel value */
2801
result.red = (MagickRealType) p[r].red;
2802
result.green = (MagickRealType) p[r].green;
2803
result.blue = (MagickRealType) p[r].blue;
2804
result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
2805
result.index = 0.0;
2806
if ( image->colorspace == CMYKColorspace)
2807
result.index = (MagickRealType) p_indexes[r];
2808
1188
result.green=(MagickRealType) 0;
1189
result.blue=(MagickRealType) 0;
1190
result.opacity=(MagickRealType) 0;
1191
result.index=(MagickRealType) 0;
2809
1192
switch (method) {
2810
1193
case ConvolveMorphology:
2811
/* Set the bias of the weighted average output */
2812
result.red =
2813
result.green =
2814
result.blue =
2815
result.opacity =
2816
result.index = bias;
1194
/* Set the user defined bias of the weighted average output
1195
**
1196
** FUTURE: provide some way for internal functions to disable
1197
** user defined bias and scaling effects.
1198
*/
1199
result=bias;
1200
break;
1201
case DilateMorphology:
1202
result.red =
1203
result.green =
1204
result.blue =
1205
result.opacity =
1206
result.index = -MagickHuge;
1207
break;
1208
case ErodeMorphology:
1209
result.red =
1210
result.green =
1211
result.blue =
1212
result.opacity =
1213
result.index = +MagickHuge;
2817
1214
break;
2818
1215
case DilateIntensityMorphology:
2819
1216
case ErodeIntensityMorphology:
2820
/* use a boolean flag indicating when first match found */
2821
result.red = 0.0; /* result is not used otherwise */
1217
result.red = 0.0; /* flag indicating first match found */
2822
1218
break;
2823
1219
default:
1220
/* Otherwise just start with the original pixel value */
1221
result.red = (MagickRealType) p[r].red;
1222
result.green = (MagickRealType) p[r].green;
1223
result.blue = (MagickRealType) p[r].blue;
1224
result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1225
if ( image->colorspace == CMYKColorspace)
1226
result.index = (MagickRealType) p_indexes[r];
2824
1227
break;
2825
1228
}
2826
1229
2836
1239
** the kernel, and thus 'lower-level' that Convolution. However
2837
1240
** as Convolution is the more common method used, and it does not
2838
1241
** really cost us much in terms of processing to use a reflected
2839
** kernel, so it is Convolution that is implemented.
1242
** kernel it is Convolution that is implemented.
2840
1243
**
2841
1244
** Correlation will have its kernel reflected before calling
2842
1245
** this function to do a Convolve.
2844
1247
** For more details of Correlation vs Convolution see
2845
1248
** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2846
1249
*/
2847
k = &kernel->values[ kernel->width*kernel->height-1 ];
2848
k_pixels = p;
2849
k_indexes = p_indexes;
2850
if ( ((channel & SyncChannels) == 0 ) ||
2851
(image->matte == MagickFalse) )
2852
{ /* No 'Sync' involved.
2853
** Convolution is simple greyscale channel operation
2854
*/
2855
for (v=0; v < (ssize_t) kernel->height; v++) {
2856
for (u=0; u < (ssize_t) kernel->width; u++, k--) {
1250
if (((channel & OpacityChannel) == 0) ||
1251
(image->matte == MagickFalse))
1252
{
1253
/* Convolution without transparency effects */
1254
k = &kernel->values[ kernel->width*kernel->height-1 ];
1255
k_pixels = p;
1256
k_indexes = p_indexes;
1257
for (v=0; v < (long) kernel->height; v++) {
1258
for (u=0; u < (long) kernel->width; u++, k--) {
2857
1259
if ( IsNan(*k) ) continue;
2858
1260
result.red += (*k)*k_pixels[u].red;
2859
1261
result.green += (*k)*k_pixels[u].green;
2860
1262
result.blue += (*k)*k_pixels[u].blue;
2861
result.opacity += (*k)*k_pixels[u].opacity;
1263
/* result.opacity += not involved here */
2862
1264
if ( image->colorspace == CMYKColorspace)
2863
1265
result.index += (*k)*k_indexes[u];
2864
1266
}
2865
1267
k_pixels += image->columns+kernel->width;
2866
1268
k_indexes += image->columns+kernel->width;
2867
1269
}
2868
if ((channel & RedChannel) != 0)
2869
q->red = ClampToQuantum(result.red);
2870
if ((channel & GreenChannel) != 0)
2871
q->green = ClampToQuantum(result.green);
2872
if ((channel & BlueChannel) != 0)
2873
q->blue = ClampToQuantum(result.blue);
2874
if ((channel & OpacityChannel) != 0
2875
&& image->matte == MagickTrue )
2876
q->opacity = ClampToQuantum(result.opacity);
2877
if ((channel & IndexChannel) != 0
2878
&& image->colorspace == CMYKColorspace)
2879
q_indexes[x] = ClampToQuantum(result.index);
2880
1270
}
2881
1271
else
2882
{ /* Channel 'Sync' Flag, and Alpha Channel enabled.
2883
** Weight the color channels with Alpha Channel so that
2884
** transparent pixels are not part of the results.
2885
*/
1272
{ /* Kernel & Alpha weighted Convolution */
2886
1273
MagickRealType
2887
alpha, /* alpha weighting of colors : kernel*alpha */
2888
gamma; /* divisor, sum of color weighting values */
1274
alpha, /* alpha value * kernel weighting */
1275
gamma; /* weighting divisor */
2889
1276
2890
1277
gamma=0.0;
2891
for (v=0; v < (ssize_t) kernel->height; v++) {
2892
for (u=0; u < (ssize_t) kernel->width; u++, k--) {
1278
k = &kernel->values[ kernel->width*kernel->height-1 ];
1279
k_pixels = p;
1280
k_indexes = p_indexes;
1281
for (v=0; v < (long) kernel->height; v++) {
1282
for (u=0; u < (long) kernel->width; u++, k--) {
2893
1283
if ( IsNan(*k) ) continue;
2894
1284
alpha=(*k)*(QuantumScale*(QuantumRange-
2895
1285
k_pixels[u].opacity));
2897
1287
result.red += alpha*k_pixels[u].red;
2898
1288
result.green += alpha*k_pixels[u].green;
2899
1289
result.blue += alpha*k_pixels[u].blue;
2900
result.opacity += (*k)*k_pixels[u].opacity;
1290
result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
2901
1291
if ( image->colorspace == CMYKColorspace)
2902
1292
result.index += alpha*k_indexes[u];
2903
1293
}
2904
1294
k_pixels += image->columns+kernel->width;
2905
1295
k_indexes += image->columns+kernel->width;
2906
1296
}
2907
/* Sync'ed channels, all channels are modified */
2908
1297
gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2909
q->red = ClampToQuantum(gamma*result.red);
2910
q->green = ClampToQuantum(gamma*result.green);
2911
q->blue = ClampToQuantum(gamma*result.blue);
2912
q->opacity = ClampToQuantum(result.opacity);
2913
if (image->colorspace == CMYKColorspace)
2914
q_indexes[x] = ClampToQuantum(gamma*result.index);
1298
result.red *= gamma;
1299
result.green *= gamma;
1300
result.blue *= gamma;
1301
result.opacity *= gamma;
1302
result.index *= gamma;
2915
1303
}
2916
1304
break;
2917
1305
2918
1306
case ErodeMorphology:
2919
/* Minimum Value within kernel neighbourhood
1307
/* Minimize Value within kernel neighbourhood
2920
1308
**
2921
1309
** NOTE that the kernel is not reflected for this operation!
2922
1310
**
2927
1315
k = kernel->values;
2928
1316
k_pixels = p;
2929
1317
k_indexes = p_indexes;
2930
for (v=0; v < (ssize_t) kernel->height; v++) {
2931
for (u=0; u < (ssize_t) kernel->width; u++, k++) {
1318
for (v=0; v < (long) kernel->height; v++) {
1319
for (u=0; u < (long) kernel->width; u++, k++) {
2932
1320
if ( IsNan(*k) || (*k) < 0.5 ) continue;
2933
Minimize(min.red, (double) k_pixels[u].red);
2934
Minimize(min.green, (double) k_pixels[u].green);
2935
Minimize(min.blue, (double) k_pixels[u].blue);
2936
Minimize(min.opacity,
2937
QuantumRange-(double) k_pixels[u].opacity);
1321
Minimize(result.red, (double) k_pixels[u].red);
1322
Minimize(result.green, (double) k_pixels[u].green);
1323
Minimize(result.blue, (double) k_pixels[u].blue);
1324
Minimize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
2938
1325
if ( image->colorspace == CMYKColorspace)
2939
Minimize(min.index, (double) k_indexes[u]);
1326
Minimize(result.index, (double) k_indexes[u]);
2940
1327
}
2941
1328
k_pixels += image->columns+kernel->width;
2942
1329
k_indexes += image->columns+kernel->width;
2944
1331
break;
2945
1332
2946
1333
case DilateMorphology:
2947
/* Maximum Value within kernel neighbourhood
1334
/* Maximize Value within kernel neighbourhood
2948
1335
**
2949
1336
** NOTE for correct working of this operation for asymetrical
2950
1337
** kernels, the kernel needs to be applied in its reflected form.
2958
1345
k = &kernel->values[ kernel->width*kernel->height-1 ];
2959
1346
k_pixels = p;
2960
1347
k_indexes = p_indexes;
2961
for (v=0; v < (ssize_t) kernel->height; v++) {
2962
for (u=0; u < (ssize_t) kernel->width; u++, k--) {
1348
for (v=0; v < (long) kernel->height; v++) {
1349
for (u=0; u < (long) kernel->width; u++, k--) {
2963
1350
if ( IsNan(*k) || (*k) < 0.5 ) continue;
2964
Maximize(max.red, (double) k_pixels[u].red);
2965
Maximize(max.green, (double) k_pixels[u].green);
2966
Maximize(max.blue, (double) k_pixels[u].blue);
2967
Maximize(max.opacity,
2968
QuantumRange-(double) k_pixels[u].opacity);
1351
Maximize(result.red, (double) k_pixels[u].red);
1352
Maximize(result.green, (double) k_pixels[u].green);
1353
Maximize(result.blue, (double) k_pixels[u].blue);
1354
Maximize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
2969
1355
if ( image->colorspace == CMYKColorspace)
2970
Maximize(max.index, (double) k_indexes[u]);
2971
}
2972
k_pixels += image->columns+kernel->width;
2973
k_indexes += image->columns+kernel->width;
2974
}
2975
break;
2976
2977
case HitAndMissMorphology:
2978
case ThinningMorphology:
2979
case ThickenMorphology:
2980
/* Minimum of Foreground Pixel minus Maxumum of Background Pixels
2981
**
2982
** NOTE that the kernel is not reflected for this operation,
2983
** and consists of both foreground and background pixel
2984
** neighbourhoods, 0.0 for background, and 1.0 for foreground
2985
** with either Nan or 0.5 values for don't care.
2986
**
2987
** Note that this will never produce a meaningless negative
2988
** result. Such results can cause Thinning/Thicken to not work
2989
** correctly when used against a greyscale image.
2990
*/
2991
k = kernel->values;
2992
k_pixels = p;
2993
k_indexes = p_indexes;
2994
for (v=0; v < (ssize_t) kernel->height; v++) {
2995
for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2996
if ( IsNan(*k) ) continue;
2997
if ( (*k) > 0.7 )
2998
{ /* minimim of foreground pixels */
2999
Minimize(min.red, (double) k_pixels[u].red);
3000
Minimize(min.green, (double) k_pixels[u].green);
3001
Minimize(min.blue, (double) k_pixels[u].blue);
3002
Minimize(min.opacity,
3003
QuantumRange-(double) k_pixels[u].opacity);
3004
if ( image->colorspace == CMYKColorspace)
3005
Minimize(min.index, (double) k_indexes[u]);
3006
}
3007
else if ( (*k) < 0.3 )
3008
{ /* maximum of background pixels */
3009
Maximize(max.red, (double) k_pixels[u].red);
3010
Maximize(max.green, (double) k_pixels[u].green);
3011
Maximize(max.blue, (double) k_pixels[u].blue);
3012
Maximize(max.opacity,
3013
QuantumRange-(double) k_pixels[u].opacity);
3014
if ( image->colorspace == CMYKColorspace)
3015
Maximize(max.index, (double) k_indexes[u]);
3016
}
3017
}
3018
k_pixels += image->columns+kernel->width;
3019
k_indexes += image->columns+kernel->width;
3020
}
3021
/* Pattern Match if difference is positive */
3022
min.red -= max.red; Maximize( min.red, 0.0 );
3023
min.green -= max.green; Maximize( min.green, 0.0 );
3024
min.blue -= max.blue; Maximize( min.blue, 0.0 );
3025
min.opacity -= max.opacity; Maximize( min.opacity, 0.0 );
3026
min.index -= max.index; Maximize( min.index, 0.0 );
1356
Maximize(result.index, (double) k_indexes[u]);
1357
}
1358
k_pixels += image->columns+kernel->width;
1359
k_indexes += image->columns+kernel->width;
1360
}
3027
1361
break;
3028
1362
3029
1363
case ErodeIntensityMorphology:
3030
1364
/* Select Pixel with Minimum Intensity within kernel neighbourhood
3031
1365
**
3032
1366
** WARNING: the intensity test fails for CMYK and does not
3033
** take into account the moderating effect of the alpha channel
1367
** take into account the moderating effect of teh alpha channel
3034
1368
** on the intensity.
3035
1369
**
3036
1370
** NOTE that the kernel is not reflected for this operation!
3038
1372
k = kernel->values;
3039
1373
k_pixels = p;
3040
1374
k_indexes = p_indexes;
3041
for (v=0; v < (ssize_t) kernel->height; v++) {
3042
for (u=0; u < (ssize_t) kernel->width; u++, k++) {
1375
for (v=0; v < (long) kernel->height; v++) {
1376
for (u=0; u < (long) kernel->width; u++, k++) {
3043
1377
if ( IsNan(*k) || (*k) < 0.5 ) continue;
3044
1378
if ( result.red == 0.0 ||
3045
1379
PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
3058
1392
/* Select Pixel with Maximum Intensity within kernel neighbourhood
3059
1393
**
3060
1394
** WARNING: the intensity test fails for CMYK and does not
3061
** take into account the moderating effect of the alpha channel
3062
** on the intensity (yet).
1395
** take into account the moderating effect of teh alpha channel
1396
** on the intensity.
3063
1397
**
3064
1398
** NOTE for correct working of this operation for asymetrical
3065
1399
** kernels, the kernel needs to be applied in its reflected form.
3068
1402
k = &kernel->values[ kernel->width*kernel->height-1 ];
3069
1403
k_pixels = p;
3070
1404
k_indexes = p_indexes;
3071
for (v=0; v < (ssize_t) kernel->height; v++) {
3072
for (u=0; u < (ssize_t) kernel->width; u++, k--) {
1405
for (v=0; v < (long) kernel->height; v++) {
1406
for (u=0; u < (long) kernel->width; u++, k--) {
3073
1407
if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3074
1408
if ( result.red == 0.0 ||
3075
1409
PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
3084
1418
}
3085
1419
break;
3086
1420
3087
3088
1421
case DistanceMorphology:
3089
1422
/* Add kernel Value and select the minimum value found.
3090
1423
** The result is a iterative distance from edge of image shape.
3093
1426
** be the case. For example how about a distance from left edges?
3094
1427
** To work correctly with asymetrical kernels the reflected kernel
3095
1428
** needs to be applied.
3096
**
3097
** Actually this is really a GreyErode with a negative kernel!
3098
**
3099
*/
1429
*/
1430
#if 0
1431
/* No need to do distance morphology if original value is zero
1432
** Unfortunatally I have not been able to get this right
1433
** when channel selection also becomes involved. -- Arrgghhh
1434
*/
1435
if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1436
|| ((channel & GreenChannel) == 0 && p[r].green == 0)
1437
|| ((channel & BlueChannel) == 0 && p[r].blue == 0)
1438
|| ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1439
|| (( (channel & IndexChannel) == 0
1440
|| image->colorspace != CMYKColorspace
1441
) && p_indexes[x] ==0 )
1442
)
1443
break;
1444
#endif
3100
1445
k = &kernel->values[ kernel->width*kernel->height-1 ];
3101
1446
k_pixels = p;
3102
1447
k_indexes = p_indexes;
3103
for (v=0; v < (ssize_t) kernel->height; v++) {
3104
for (u=0; u < (ssize_t) kernel->width; u++, k--) {
1448
for (v=0; v < (long) kernel->height; v++) {
1449
for (u=0; u < (long) kernel->width; u++, k--) {
3105
1450
if ( IsNan(*k) ) continue;
3106
1451
Minimize(result.red, (*k)+k_pixels[u].red);
3107
1452
Minimize(result.green, (*k)+k_pixels[u].green);
3119
1464
default:
3120
1465
break; /* Do nothing */
3121
1466
}
3122
/* Final mathematics of results (combine with original image?)
3123
**
3124
** NOTE: Difference Morphology operators Edge* and *Hat could also
3125
** be done here but works better with iteration as a image difference
3126
** in the controling function (below). Thicken and Thinning however
3127
** should be done here so thay can be iterated correctly.
3128
*/
3129
switch ( method ) {
3130
case HitAndMissMorphology:
3131
case ErodeMorphology:
3132
result = min; /* minimum of neighbourhood */
3133
break;
3134
case DilateMorphology:
3135
result = max; /* maximum of neighbourhood */
3136
break;
3137
case ThinningMorphology:
3138
/* subtract pattern match from original */
3139
result.red -= min.red;
3140
result.green -= min.green;
3141
result.blue -= min.blue;
3142
result.opacity -= min.opacity;
3143
result.index -= min.index;
3144
break;
3145
case ThickenMorphology:
3146
/* Add the pattern matchs to the original */
3147
result.red += min.red;
3148
result.green += min.green;
3149
result.blue += min.blue;
3150
result.opacity += min.opacity;
3151
result.index += min.index;
3152
break;
3153
default:
3154
/* result directly calculated or assigned */
3155
break;
3156
}
3157
/* Assign the resulting pixel values - Clamping Result */
3158
1467
switch ( method ) {
3159
1468
case UndefinedMorphology:
3160
case ConvolveMorphology:
3161
1469
case DilateIntensityMorphology:
3162
1470
case ErodeIntensityMorphology:
3163
1471
break; /* full pixel was directly assigned - not a channel method */
3164
1472
default:
1473
/* Assign the results */
3165
1474
if ((channel & RedChannel) != 0)
3166
1475
q->red = ClampToQuantum(result.red);
3167
1476
if ((channel & GreenChannel) != 0)
3176
1485
q_indexes[x] = ClampToQuantum(result.index);
3177
1486
break;
3178
1487
}
3179
/* Count up changed pixels */
3180
1488
if ( ( p[r].red != q->red )
3181
1489
|| ( p[r].green != q->green )
3182
1490
|| ( p[r].blue != q->blue )
3183
1491
|| ( p[r].opacity != q->opacity )
3184
1492
|| ( image->colorspace == CMYKColorspace &&
3185
1493
p_indexes[r] != q_indexes[x] ) )
3186
changed++; /* The pixel was changed in some way! */
1494
changed++; /* The pixel had some value changed! */
3187
1495
p++;
3188
1496
q++;
3189
1497
} /* x */
3190
if ( SyncCacheViewAuthenticPixels(q_view,exception) == MagickFalse)
1498
sync=SyncCacheViewAuthenticPixels(q_view,exception);
1499
if (sync == MagickFalse)
3191
1500
status=MagickFalse;
3192
1501
if (image->progress_monitor != (MagickProgressMonitor) NULL)
3193
1502
{
3205
1514
result_image->type=image->type;
3206
1515
q_view=DestroyCacheView(q_view);
3207
1516
p_view=DestroyCacheView(p_view);
3208
return(status ? (size_t) changed : 0);
3209
}
3210
3211
3212
MagickExport Image *MorphologyApply(const Image *image, const ChannelType
3213
channel,const MorphologyMethod method, const ssize_t iterations,
3214
const KernelInfo *kernel, const CompositeOperator compose,
3215
const double bias, ExceptionInfo *exception)
3216
{
3217
Image
3218
*curr_image, /* Image we are working with or iterating */
3219
*work_image, /* secondary image for primative iteration */
3220
*save_image, /* saved image - for 'edge' method only */
3221
*rslt_image; /* resultant image - after multi-kernel handling */
1517
return(status ? (unsigned long) changed : 0);
1518
}
1519
1520
1521
MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1522
method, const long iterations,const KernelInfo *kernel, ExceptionInfo
1523
*exception)
1524
{
1525
Image
1526
*morphology_image;
1527
1528
morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1529
iterations,kernel,exception);
1530
return(morphology_image);
1531
}
1532
1533
1534
MagickExport Image *MorphologyImageChannel(const Image *image,
1535
const ChannelType channel,const MorphologyMethod method,
1536
const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
1537
{
1538
long
1539
count;
1540
1541
Image
1542
*new_image,
1543
*old_image,
1544
*grad_image;
1545
1546
const char
1547
*artifact;
1548
1549
unsigned long
1550
changed,
1551
limit;
3222
1552
3223
1553
KernelInfo
3224
*reflected_kernel, /* A reflected copy of the kernel (if needed) */
3225
*norm_kernel, /* the current normal un-reflected kernel */
3226
*rflt_kernel, /* the current reflected kernel (if needed) */
3227
*this_kernel; /* the kernel being applied */
1554
*curr_kernel;
3228
1555
3229
1556
MorphologyMethod
3230
primative; /* the current morphology primative being applied */
3231
3232
CompositeOperator
3233
rslt_compose; /* multi-kernel compose method for results to use */
3234
3235
MagickBooleanType
3236
verbose; /* verbose output of results */
3237
3238
size_t
3239
method_loop, /* Loop 1: number of compound method iterations */
3240
method_limit, /* maximum number of compound method iterations */
3241
kernel_number, /* Loop 2: the kernel number being applied */
3242
stage_loop, /* Loop 3: primative loop for compound morphology */
3243
stage_limit, /* how many primatives in this compound */
3244
kernel_loop, /* Loop 4: iterate the kernel (basic morphology) */
3245
kernel_limit, /* number of times to iterate kernel */
3246
count, /* total count of primative steps applied */
3247
changed, /* number pixels changed by last primative operation */
3248
kernel_changed, /* total count of changed using iterated kernel */
3249
method_changed; /* total count of changed over method iteration */
3250
3251
char
3252
v_info[80];
1557
curr_method;
3253
1558
3254
1559
assert(image != (Image *) NULL);
3255
1560
assert(image->signature == MagickSignature);
3258
1563
assert(exception != (ExceptionInfo *) NULL);
3259
1564
assert(exception->signature == MagickSignature);
3260
1565
3261
count = 0; /* number of low-level morphology primatives performed */
3262
1566
if ( iterations == 0 )
3263
return((Image *)NULL); /* null operation - nothing to do! */
3264
3265
kernel_limit = (size_t) iterations;
3266
if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3267
kernel_limit = image->columns > image->rows ? image->columns : image->rows;
3268
3269
verbose = ( GetImageArtifact(image,"verbose") != (const char *) NULL ) ?
3270
MagickTrue : MagickFalse;
3271
3272
/* initialise for cleanup */
3273
curr_image = (Image *) image;
3274
work_image = save_image = rslt_image = (Image *) NULL;
3275
reflected_kernel = (KernelInfo *) NULL;
3276
3277
/* Initialize specific methods
3278
* + which loop should use the given iteratations
3279
* + how many primatives make up the compound morphology
3280
* + multi-kernel compose method to use (by default)
1567
return((Image *)NULL); /* null operation - nothing to do! */
1568
1569
/* kernel must be valid at this point
1570
* (except maybe for posible future morphology methods like "Prune"
3281
1571
*/
3282
method_limit = 1; /* just do method once, unless otherwise set */
3283
stage_limit = 1; /* assume method is not a compount */
3284
rslt_compose = compose; /* and we are composing multi-kernels as given */
3285
switch( method ) {
3286
case SmoothMorphology: /* 4 primative compound morphology */
3287
stage_limit = 4;
3288
break;
3289
case OpenMorphology: /* 2 primative compound morphology */
1572
assert(kernel != (KernelInfo *)NULL);
1573
1574
count = 0; /* interation count */
1575
changed = 1; /* if compound method assume image was changed */
1576
curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
1577
curr_method = method; /* to be changed as nessary */
1578
1579
limit = (unsigned long) iterations;
1580
if ( iterations < 0 )
1581
limit = image->columns > image->rows ? image->columns : image->rows;
1582
1583
/* Third-level morphology methods */
1584
grad_image=(Image *) NULL;
1585
switch( curr_method ) {
1586
case EdgeMorphology:
1587
grad_image = MorphologyImageChannel(image, channel,
1588
DilateMorphology, iterations, curr_kernel, exception);
1589
/* FALL-THRU */
1590
case EdgeInMorphology:
1591
curr_method = ErodeMorphology;
1592
break;
1593
case EdgeOutMorphology:
1594
curr_method = DilateMorphology;
1595
break;
1596
case TopHatMorphology:
1597
curr_method = OpenMorphology;
1598
break;
1599
case BottomHatMorphology:
1600
curr_method = CloseMorphology;
1601
break;
1602
default:
1603
break; /* not a third-level method */
1604
}
1605
1606
/* Second-level morphology methods */
1607
switch( curr_method ) {
1608
case OpenMorphology:
1609
/* Open is a Erode then a Dilate without reflection */
1610
new_image = MorphologyImageChannel(image, channel,
1611
ErodeMorphology, iterations, curr_kernel, exception);
1612
if (new_image == (Image *) NULL)
1613
return((Image *) NULL);
1614
curr_method = DilateMorphology;
1615
break;
3290
1616
case OpenIntensityMorphology:
3291
case TopHatMorphology:
1617
new_image = MorphologyImageChannel(image, channel,
1618
ErodeIntensityMorphology, iterations, curr_kernel, exception);
1619
if (new_image == (Image *) NULL)
1620
return((Image *) NULL);
1621
curr_method = DilateIntensityMorphology;
1622
break;
1623
3292
1624
case CloseMorphology:
1625
/* Close is a Dilate then Erode using reflected kernel */
1626
/* A reflected kernel is needed for a Close */
1627
if ( curr_kernel == kernel )
1628
curr_kernel = CloneKernelInfo(kernel);
1629
RotateKernelInfo(curr_kernel,180);
1630
new_image = MorphologyImageChannel(image, channel,
1631
DilateMorphology, iterations, curr_kernel, exception);
1632
if (new_image == (Image *) NULL)
1633
return((Image *) NULL);
1634
curr_method = ErodeMorphology;
1635
break;
3293
1636
case CloseIntensityMorphology:
3294
case BottomHatMorphology:
3295
case EdgeMorphology:
3296
stage_limit = 2;
3297
break;
3298
case HitAndMissMorphology:
3299
rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3300
/* FALL THUR */
3301
case ThinningMorphology:
3302
case ThickenMorphology:
3303
method_limit = kernel_limit; /* iterate the whole method */
3304
kernel_limit = 1; /* do not do kernel iteration */
3305
break;
3306
default:
3307
break;
3308
}
3309
3310
/* Handle user (caller) specified multi-kernel composition method */
3311
if ( compose != UndefinedCompositeOp )
3312
rslt_compose = compose; /* override default composition for method */
3313
if ( rslt_compose == UndefinedCompositeOp )
3314
rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3315
3316
/* Some methods require a reflected kernel to use with primatives.
3317
* Create the reflected kernel for those methods. */
3318
switch ( method ) {
1637
/* A reflected kernel is needed for a Close */
1638
if ( curr_kernel == kernel )
1639
curr_kernel = CloneKernelInfo(kernel);
1640
RotateKernelInfo(curr_kernel,180);
1641
new_image = MorphologyImageChannel(image, channel,
1642
DilateIntensityMorphology, iterations, curr_kernel, exception);
1643
if (new_image == (Image *) NULL)
1644
return((Image *) NULL);
1645
curr_method = ErodeIntensityMorphology;
1646
break;
1647
3319
1648
case CorrelateMorphology:
3320
case CloseMorphology:
3321
case CloseIntensityMorphology:
3322
case BottomHatMorphology:
3323
case SmoothMorphology:
3324
reflected_kernel = CloneKernelInfo(kernel);
3325
if (reflected_kernel == (KernelInfo *) NULL)
3326
goto error_cleanup;
3327
RotateKernelInfo(reflected_kernel,180);
3328
break;
3329
default:
3330
break;
3331
}
3332
3333
/* Loop 1: iterate the compound method */
3334
method_loop = 0;
3335
method_changed = 1;
3336
while ( method_loop < method_limit && method_changed > 0 ) {
3337
method_loop++;
3338
method_changed = 0;
3339
3340
/* Loop 2: iterate over each kernel in a multi-kernel list */
3341
norm_kernel = (KernelInfo *) kernel;
3342
this_kernel = (KernelInfo *) kernel;
3343
rflt_kernel = reflected_kernel;
3344
3345
kernel_number = 0;
3346
while ( norm_kernel != NULL ) {
3347
3348
/* Loop 3: Compound Morphology Staging - Select Primative to apply */
3349
stage_loop = 0; /* the compound morphology stage number */
3350
while ( stage_loop < stage_limit ) {
3351
stage_loop++; /* The stage of the compound morphology */
3352
3353
/* Select primative morphology for this stage of compound method */
3354
this_kernel = norm_kernel; /* default use unreflected kernel */
3355
primative = method; /* Assume method is a primative */
3356
switch( method ) {
3357
case ErodeMorphology: /* just erode */
3358
case EdgeInMorphology: /* erode and image difference */
3359
primative = ErodeMorphology;
3360
break;
3361
case DilateMorphology: /* just dilate */
3362
case EdgeOutMorphology: /* dilate and image difference */
3363
primative = DilateMorphology;
3364
break;
3365
case OpenMorphology: /* erode then dialate */
3366
case TopHatMorphology: /* open and image difference */
3367
primative = ErodeMorphology;
3368
if ( stage_loop == 2 )
3369
primative = DilateMorphology;
3370
break;
3371
case OpenIntensityMorphology:
3372
primative = ErodeIntensityMorphology;
3373
if ( stage_loop == 2 )
3374
primative = DilateIntensityMorphology;
3375
break;
3376
case CloseMorphology: /* dilate, then erode */
3377
case BottomHatMorphology: /* close and image difference */
3378
this_kernel = rflt_kernel; /* use the reflected kernel */
3379
primative = DilateMorphology;
3380
if ( stage_loop == 2 )
3381
primative = ErodeMorphology;
3382
break;
3383
case CloseIntensityMorphology:
3384
this_kernel = rflt_kernel; /* use the reflected kernel */
3385
primative = DilateIntensityMorphology;
3386
if ( stage_loop == 2 )
3387
primative = ErodeIntensityMorphology;
3388
break;
3389
case SmoothMorphology: /* open, close */
3390
switch ( stage_loop ) {
3391
case 1: /* start an open method, which starts with Erode */
3392
primative = ErodeMorphology;
3393
break;
3394
case 2: /* now Dilate the Erode */
3395
primative = DilateMorphology;
3396
break;
3397
case 3: /* Reflect kernel a close */
3398
this_kernel = rflt_kernel; /* use the reflected kernel */
3399
primative = DilateMorphology;
3400
break;
3401
case 4: /* Finish the Close */
3402
this_kernel = rflt_kernel; /* use the reflected kernel */
3403
primative = ErodeMorphology;
3404
break;
3405
}
3406
break;
3407
case EdgeMorphology: /* dilate and erode difference */
3408
primative = DilateMorphology;
3409
if ( stage_loop == 2 ) {
3410
save_image = curr_image; /* save the image difference */
3411
curr_image = (Image *) image;
3412
primative = ErodeMorphology;
3413
}
3414
break;
3415
case CorrelateMorphology:
3416
/* A Correlation is a Convolution with a reflected kernel.
3417
** However a Convolution is a weighted sum using a reflected
3418
** kernel. It may seem stange to convert a Correlation into a
3419
** Convolution as the Correlation is the simplier method, but
3420
** Convolution is much more commonly used, and it makes sense to
3421
** implement it directly so as to avoid the need to duplicate the
3422
** kernel when it is not required (which is typically the
3423
** default).
3424
*/
3425
this_kernel = rflt_kernel; /* use the reflected kernel */
3426
primative = ConvolveMorphology;
3427
break;
3428
default:
3429
break;
3430
}
3431
assert( this_kernel != (KernelInfo *) NULL );
3432
3433
/* Extra information for debugging compound operations */
3434
if ( verbose == MagickTrue ) {
3435
if ( stage_limit > 1 )
3436
(void) FormatMagickString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
3437
MagickOptionToMnemonic(MagickMorphologyOptions,method),(double)
3438
method_loop,(double) stage_loop);
3439
else if ( primative != method )
3440
(void) FormatMagickString(v_info, MaxTextExtent, "%s:%.20g -> ",
3441
MagickOptionToMnemonic(MagickMorphologyOptions, method),(double)
3442
method_loop);
3443
else
3444
v_info[0] = '\0';
3445
}
3446
3447
/* Loop 4: Iterate the kernel with primative */
3448
kernel_loop = 0;
3449
kernel_changed = 0;
3450
changed = 1;
3451
while ( kernel_loop < kernel_limit && changed > 0 ) {
3452
kernel_loop++; /* the iteration of this kernel */
3453
3454
/* Create a destination image, if not yet defined */
3455
if ( work_image == (Image *) NULL )
3456
{
3457
work_image=CloneImage(image,0,0,MagickTrue,exception);
3458
if (work_image == (Image *) NULL)
3459
goto error_cleanup;
3460
if (SetImageStorageClass(work_image,DirectClass) == MagickFalse)
3461
{
3462
InheritException(exception,&work_image->exception);
3463
goto error_cleanup;
3464
}
3465
}
3466
3467
/* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3468
count++;
3469
changed = MorphologyPrimitive(curr_image, work_image, primative,
3470
channel, this_kernel, bias, exception);
3471
kernel_changed += changed;
3472
method_changed += changed;
3473
3474
if ( verbose == MagickTrue ) {
3475
if ( kernel_loop > 1 )
3476
fprintf(stderr, "\n"); /* add end-of-line from previous */
3477
(void) fprintf(stderr, "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3478
v_info,MagickOptionToMnemonic(MagickMorphologyOptions,
3479
primative),(this_kernel == rflt_kernel ) ? "*" : "",
3480
(double) (method_loop+kernel_loop-1),(double) kernel_number,
3481
(double) count,(double) changed);
3482
}
3483
/* prepare next loop */
3484
{ Image *tmp = work_image; /* swap images for iteration */
3485
work_image = curr_image;
3486
curr_image = tmp;
3487
}
3488
if ( work_image == image )
3489
work_image = (Image *) NULL; /* replace input 'image' */
3490
3491
} /* End Loop 4: Iterate the kernel with primative */
3492
3493
if ( verbose == MagickTrue && kernel_changed != changed )
3494
fprintf(stderr, " Total %.20g",(double) kernel_changed);
3495
if ( verbose == MagickTrue && stage_loop < stage_limit )
3496
fprintf(stderr, "\n"); /* add end-of-line before looping */
3497
3498
#if 0
3499
fprintf(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3500
fprintf(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
3501
fprintf(stderr, " work =0x%lx\n", (unsigned long)work_image);
3502
fprintf(stderr, " save =0x%lx\n", (unsigned long)save_image);
3503
fprintf(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
3504
#endif
3505
3506
} /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3507
3508
/* Final Post-processing for some Compound Methods
1649
/* A Correlation is actually a Convolution with a reflected kernel.
1650
** However a Convolution is a weighted sum with a reflected kernel.
1651
** It may seem stange to convert a Correlation into a Convolution
1652
** as the Correleation is the simplier method, but Convolution is
1653
** much more commonly used, and it makes sense to implement it directly
1654
** so as to avoid the need to duplicate the kernel when it is not
1655
** required (which is typically the default).
1656
*/
1657
if ( curr_kernel == kernel )
1658
curr_kernel = CloneKernelInfo(kernel);
1659
RotateKernelInfo(curr_kernel,180);
1660
curr_method = ConvolveMorphology;
1661
/* FALL-THRU into Correlate (weigthed sum without reflection) */
1662
1663
case ConvolveMorphology:
1664
/* Scale or Normalize kernel, according to user wishes
1665
** before using it for the Convolve/Correlate method.
3509
1666
**
3510
** The removal of any 'Sync' channel flag in the Image Compositon
3511
** below ensures the methematical compose method is applied in a
3512
** purely mathematical way, and only to the selected channels.
3513
** Turn off SVG composition 'alpha blending'.
1667
** FUTURE: provide some way for internal functions to disable
1668
** user bias and scaling effects.
3514
1669
*/
3515
switch( method ) {
3516
case EdgeOutMorphology:
3517
case EdgeInMorphology:
3518
case TopHatMorphology:
3519
case BottomHatMorphology:
3520
if ( verbose == MagickTrue )
3521
fprintf(stderr, "\n%s: Difference with original image",
3522
MagickOptionToMnemonic(MagickMorphologyOptions, method) );
3523
(void) CompositeImageChannel(curr_image,
3524
(ChannelType) (channel & ~SyncChannels),
3525
DifferenceCompositeOp, image, 0, 0);
3526
break;
3527
case EdgeMorphology:
3528
if ( verbose == MagickTrue )
3529
fprintf(stderr, "\n%s: Difference of Dilate and Erode",
3530
MagickOptionToMnemonic(MagickMorphologyOptions, method) );
3531
(void) CompositeImageChannel(curr_image,
3532
(ChannelType) (channel & ~SyncChannels),
3533
DifferenceCompositeOp, save_image, 0, 0);
3534
save_image = DestroyImage(save_image); /* finished with save image */
3535
break;
3536
default:
3537
break;
3538
}
3539
3540
/* multi-kernel handling: re-iterate, or compose results */
3541
if ( kernel->next == (KernelInfo *) NULL )
3542
rslt_image = curr_image; /* just return the resulting image */
3543
else if ( rslt_compose == NoCompositeOp )
3544
{ if ( verbose == MagickTrue ) {
3545
if ( this_kernel->next != (KernelInfo *) NULL )
3546
fprintf(stderr, " (re-iterate)");
3547
else
3548
fprintf(stderr, " (done)");
3549
}
3550
rslt_image = curr_image; /* return result, and re-iterate */
3551
}
3552
else if ( rslt_image == (Image *) NULL)
3553
{ if ( verbose == MagickTrue )
3554
fprintf(stderr, " (save for compose)");
3555
rslt_image = curr_image;
3556
curr_image = (Image *) image; /* continue with original image */
3557
}
3558
else
3559
{ /* add the new 'current' result to the composition
3560
**
3561
** The removal of any 'Sync' channel flag in the Image Compositon
3562
** below ensures the methematical compose method is applied in a
3563
** purely mathematical way, and only to the selected channels.
3564
** Turn off SVG composition 'alpha blending'.
3565
**
3566
** The compose image order is specifically so that the new image can
3567
** be subtarcted 'Minus' from the collected result, to allow you to
3568
** convert a HitAndMiss methd into a Thinning method.
3569
*/
3570
if ( verbose == MagickTrue )
3571
fprintf(stderr, " (compose \"%s\")",
3572
MagickOptionToMnemonic(MagickComposeOptions, rslt_compose) );
3573
(void) CompositeImageChannel(curr_image,
3574
(ChannelType) (channel & ~SyncChannels), rslt_compose,
3575
rslt_image, 0, 0);
3576
rslt_image = DestroyImage(rslt_image);
3577
rslt_image = curr_image;
3578
curr_image = (Image *) image; /* continue with original image */
3579
}
3580
if ( verbose == MagickTrue )
3581
fprintf(stderr, "\n");
3582
3583
/* loop to the next kernel in a multi-kernel list */
3584
norm_kernel = norm_kernel->next;
3585
if ( rflt_kernel != (KernelInfo *) NULL )
3586
rflt_kernel = rflt_kernel->next;
3587
kernel_number++;
3588
} /* End Loop 2: Loop over each kernel */
3589
3590
} /* End Loop 1: compound method interation */
3591
3592
goto exit_cleanup;
3593
3594
/* Yes goto's are bad, but it makes cleanup lot more efficient */
3595
error_cleanup:
3596
if ( curr_image != (Image *) NULL &&
3597
curr_image != rslt_image &&
3598
curr_image != image )
3599
curr_image = DestroyImage(curr_image);
3600
if ( rslt_image != (Image *) NULL )
3601
rslt_image = DestroyImage(rslt_image);
3602
exit_cleanup:
3603
if ( curr_image != (Image *) NULL &&
3604
curr_image != rslt_image &&
3605
curr_image != image )
3606
curr_image = DestroyImage(curr_image);
3607
if ( work_image != (Image *) NULL )
3608
work_image = DestroyImage(work_image);
3609
if ( save_image != (Image *) NULL )
3610
save_image = DestroyImage(save_image);
3611
if ( reflected_kernel != (KernelInfo *) NULL )
3612
reflected_kernel = DestroyKernelInfo(reflected_kernel);
3613
return(rslt_image);
3614
}
3615
3616
/*
3617
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3618
% %
3619
% %
3620
% %
3621
% M o r p h o l o g y I m a g e C h a n n e l %
3622
% %
3623
% %
3624
% %
3625
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3626
%
3627
% MorphologyImageChannel() applies a user supplied kernel to the image
3628
% according to the given mophology method.
3629
%
3630
% This function applies any and all user defined settings before calling
3631
% the above internal function MorphologyApply().
3632
%
3633
% User defined settings include...
3634
% * Output Bias for Convolution and correlation ("-bias")
3635
% * Kernel Scale/normalize settings ("-set 'option:convolve:scale'")
3636
% This can also includes the addition of a scaled unity kernel.
3637
% * Show Kernel being applied ("-set option:showkernel 1")
3638
%
3639
% The format of the MorphologyImage method is:
3640
%
3641
% Image *MorphologyImage(const Image *image,MorphologyMethod method,
3642
% const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
3643
%
3644
% Image *MorphologyImageChannel(const Image *image, const ChannelType
3645
% channel,MorphologyMethod method,const ssize_t iterations,
3646
% KernelInfo *kernel,ExceptionInfo *exception)
3647
%
3648
% A description of each parameter follows:
3649
%
3650
% o image: the image.
3651
%
3652
% o method: the morphology method to be applied.
3653
%
3654
% o iterations: apply the operation this many times (or no change).
3655
% A value of -1 means loop until no change found.
3656
% How this is applied may depend on the morphology method.
3657
% Typically this is a value of 1.
3658
%
3659
% o channel: the channel type.
3660
%
3661
% o kernel: An array of double representing the morphology kernel.
3662
% Warning: kernel may be normalized for the Convolve method.
3663
%
3664
% o exception: return any errors or warnings in this structure.
3665
%
3666
*/
3667
3668
MagickExport Image *MorphologyImageChannel(const Image *image,
3669
const ChannelType channel,const MorphologyMethod method,
3670
const ssize_t iterations,const KernelInfo *kernel,ExceptionInfo *exception)
3671
{
3672
const char
3673
*artifact;
3674
3675
KernelInfo
3676
*curr_kernel;
3677
3678
CompositeOperator
3679
compose;
3680
3681
Image
3682
*morphology_image;
3683
3684
3685
/* Apply Convolve/Correlate Normalization and Scaling Factors.
3686
* This is done BEFORE the ShowKernelInfo() function is called so that
3687
* users can see the results of the 'option:convolve:scale' option.
3688
*/
3689
curr_kernel = (KernelInfo *) kernel;
3690
if ( method == ConvolveMorphology || method == CorrelateMorphology )
3691
{
3692
1670
artifact = GetImageArtifact(image,"convolve:scale");
3693
if ( artifact != (const char *)NULL ) {
1671
if ( artifact != (char *)NULL ) {
1672
MagickStatusType
1673
flags;
1674
GeometryInfo
1675
args;
1676
3694
1677
if ( curr_kernel == kernel )
3695
1678
curr_kernel = CloneKernelInfo(kernel);
3696
if (curr_kernel == (KernelInfo *) NULL) {
3697
curr_kernel=DestroyKernelInfo(curr_kernel);
1679
1680
args.rho = 1.0;
1681
flags = ParseGeometry(artifact, &args);
1682
ScaleKernelInfo(curr_kernel, args.rho, flags);
1683
}
1684
/* FALL-THRU to do the first, and typically the only iteration */
1685
1686
default:
1687
/* Do a single iteration using the Low-Level Morphology method!
1688
** This ensures a "new_image" has been generated, but allows us to skip
1689
** the creation of 'old_image' if no more iterations are needed.
1690
**
1691
** The "curr_method" should also be set to a low-level method that is
1692
** understood by the MorphologyApply() internal function.
1693
*/
1694
new_image=CloneImage(image,0,0,MagickTrue,exception);
1695
if (new_image == (Image *) NULL)
1696
return((Image *) NULL);
1697
if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1698
{
1699
InheritException(exception,&new_image->exception);
1700
new_image=DestroyImage(new_image);
3698
1701
return((Image *) NULL);
3699
1702
}
3700
ScaleGeometryKernelInfo(curr_kernel, artifact);
3701
}
3702
}
3703
3704
/* display the (normalized) kernel via stderr */
3705
artifact = GetImageArtifact(image,"showkernel");
3706
if ( artifact == (const char *) NULL)
3707
artifact = GetImageArtifact(image,"convolve:showkernel");
3708
if ( artifact == (const char *) NULL)
3709
artifact = GetImageArtifact(image,"morphology:showkernel");
3710
if ( artifact != (const char *) NULL)
3711
ShowKernelInfo(curr_kernel);
3712
3713
/* Override the default handling of multi-kernel morphology results
3714
* If 'Undefined' use the default method
3715
* If 'None' (default for 'Convolve') re-iterate previous result
3716
* Otherwise merge resulting images using compose method given.
3717
* Default for 'HitAndMiss' is 'Lighten'.
3718
*/
3719
compose = UndefinedCompositeOp; /* use default for method */
3720
artifact = GetImageArtifact(image,"morphology:compose");
3721
if ( artifact != (const char *) NULL)
3722
compose = (CompositeOperator) ParseMagickOption(
3723
MagickComposeOptions,MagickFalse,artifact);
3724
3725
/* Apply the Morphology */
3726
morphology_image = MorphologyApply(image, channel, method, iterations,
3727
curr_kernel, compose, image->bias, exception);
3728
3729
/* Cleanup and Exit */
1703
changed = MorphologyApply(image,new_image,curr_method,channel,curr_kernel,
1704
exception);
1705
count++;
1706
if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1707
fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1708
MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1709
count, changed);
1710
break;
1711
}
1712
1713
/* At this point the "curr_method" should not only be set to a low-level
1714
** method that is understood by the MorphologyApply() internal function,
1715
** but "new_image" should now be defined, as the image to apply the
1716
** "curr_method" to.
1717
*/
1718
1719
/* Repeat the low-level morphology until count or no change reached */
1720
if ( count < (long) limit && changed > 0 ) {
1721
old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1722
if (old_image == (Image *) NULL)
1723
return(DestroyImage(new_image));
1724
if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1725
{
1726
InheritException(exception,&old_image->exception);
1727
old_image=DestroyImage(old_image);
1728
return(DestroyImage(new_image));
1729
}
1730
while( count < (long) limit && changed != 0 )
1731
{
1732
Image *tmp = old_image;
1733
old_image = new_image;
1734
new_image = tmp;
1735
changed = MorphologyApply(old_image,new_image,curr_method,channel,
1736
curr_kernel, exception);
1737
count++;
1738
if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1739
fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1740
MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1741
count, changed);
1742
}
1743
old_image=DestroyImage(old_image);
1744
}
1745
1746
/* We are finished with kernel - destroy it if we made a clone */
3730
1747
if ( curr_kernel != kernel )
3731
1748
curr_kernel=DestroyKernelInfo(curr_kernel);
3732
return(morphology_image);
3733
}
3734
3735
MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
3736
method, const ssize_t iterations,const KernelInfo *kernel, ExceptionInfo
3737
*exception)
3738
{
3739
Image
3740
*morphology_image;
3741
3742
morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
3743
iterations,kernel,exception);
3744
return(morphology_image);
1749
1750
/* Third-level Subtractive methods post-processing */
1751
switch( method ) {
1752
case EdgeOutMorphology:
1753
case EdgeInMorphology:
1754
case TopHatMorphology:
1755
case BottomHatMorphology:
1756
/* Get Difference relative to the original image */
1757
(void) CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1758
image, 0, 0);
1759
break;
1760
case EdgeMorphology: /* subtract the Erode from a Dilate */
1761
(void) CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1762
grad_image, 0, 0);
1763
grad_image=DestroyImage(grad_image);
1764
break;
1765
default:
1766
break;
1767
}
1768
1769
return(new_image);
3745
1770
}
3746
1771
3747
1772
/*
3755
1780
% %
3756
1781
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3757
1782
%
3758
% RotateKernelInfo() rotates the kernel by the angle given.
3759
%
3760
% Currently it is restricted to 90 degree angles, of either 1D kernels
3761
% or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
3762
% It will ignore usless rotations for specific 'named' built-in kernels.
1783
% RotateKernelInfo() rotates the kernel by the angle given. Currently it is
1784
% restricted to 90 degree angles, but this may be improved in the future.
3763
1785
%
3764
1786
% The format of the RotateKernelInfo method is:
3765
1787
%
3771
1793
%
3772
1794
% o angle: angle to rotate in degrees
3773
1795
%
3774
% This function is currently internal to this module only, but can be exported
3775
% to other modules if needed.
1796
% This function is only internel to this module, as it is not finalized,
1797
% especially with regard to non-orthogonal angles, and rotation of larger
1798
% 2D kernels.
3776
1799
*/
3777
1800
static void RotateKernelInfo(KernelInfo *kernel, double angle)
3778
1801
{
3779
/* angle the lower kernels first */
3780
if ( kernel->next != (KernelInfo *) NULL)
3781
RotateKernelInfo(kernel->next, angle);
3782
3783
1802
/* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
3784
1803
**
3785
1804
** TODO: expand beyond simple 90 degree rotates, flips and flops
3790
1809
if ( angle < 0 )
3791
1810
angle += 360.0;
3792
1811
3793
if ( 337.5 < angle || angle <= 22.5 )
3794
return; /* Near zero angle - no change! - At least not at this time */
1812
if ( 315.0 < angle || angle <= 45.0 )
1813
return; /* no change! - At least at this time */
3795
1814
3796
/* Handle special cases */
3797
1815
switch (kernel->type) {
3798
1816
/* These built-in kernels are cylindrical kernels, rotating is useless */
3799
1817
case GaussianKernel:
3800
case DoGKernel:
3801
case LoGKernel:
1818
case LaplacianKernel:
1819
case LOGKernel:
1820
case DOGKernel:
3802
1821
case DiskKernel:
3803
case PeaksKernel:
3804
case LaplacianKernel:
3805
1822
case ChebyshevKernel:
3806
case ManhattanKernel:
1823
case ManhattenKernel:
3807
1824
case EuclideanKernel:
3808
1825
return;
3809
1826
3812
1829
case SquareKernel:
3813
1830
case DiamondKernel:
3814
1831
case PlusKernel:
3815
case CrossKernel:
3816
1832
return;
3817
1833
3818
1834
/* These only allows a +/-90 degree rotation (by transpose) */
3825
1841
angle -= 180;
3826
1842
break;
3827
1843
3828
default:
1844
/* these are freely rotatable in 90 degree units */
1845
case CometKernel:
1846
case UndefinedKernel:
1847
case UserDefinedKernel:
3829
1848
break;
3830
1849
}
3831
/* Attempt rotations by 45 degrees */
3832
if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
3833
{
3834
if ( kernel->width == 3 && kernel->height == 3 )
3835
{ /* Rotate a 3x3 square by 45 degree angle */
3836
MagickRealType t = kernel->values[0];
3837
kernel->values[0] = kernel->values[3];
3838
kernel->values[3] = kernel->values[6];
3839
kernel->values[6] = kernel->values[7];
3840
kernel->values[7] = kernel->values[8];
3841
kernel->values[8] = kernel->values[5];
3842
kernel->values[5] = kernel->values[2];
3843
kernel->values[2] = kernel->values[1];
3844
kernel->values[1] = t;
3845
/* rotate non-centered origin */
3846
if ( kernel->x != 1 || kernel->y != 1 ) {
3847
ssize_t x,y;
3848
x = (ssize_t) kernel->x-1;
3849
y = (ssize_t) kernel->y-1;
3850
if ( x == y ) x = 0;
3851
else if ( x == 0 ) x = -y;
3852
else if ( x == -y ) y = 0;
3853
else if ( y == 0 ) y = x;
3854
kernel->x = (ssize_t) x+1;
3855
kernel->y = (ssize_t) y+1;
3856
}
3857
angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
3858
kernel->angle = fmod(kernel->angle+45.0, 360.0);
3859
}
3860
else
3861
perror("Unable to rotate non-3x3 kernel by 45 degrees");
3862
}
3863
if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
3864
{
3865
if ( kernel->width == 1 || kernel->height == 1 )
3866
{ /* Do a transpose of a 1 dimentional kernel,
3867
** which results in a fast 90 degree rotation of some type.
3868
*/
3869
ssize_t
3870
t;
3871
t = (ssize_t) kernel->width;
3872
kernel->width = kernel->height;
3873
kernel->height = (size_t) t;
3874
t = kernel->x;
3875
kernel->x = kernel->y;
3876
kernel->y = t;
3877
if ( kernel->width == 1 ) {
3878
angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
3879
kernel->angle = fmod(kernel->angle+90.0, 360.0);
3880
} else {
3881
angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
3882
kernel->angle = fmod(kernel->angle+270.0, 360.0);
3883
}
3884
}
3885
else if ( kernel->width == kernel->height )
3886
{ /* Rotate a square array of values by 90 degrees */
3887
{ register size_t
3888
i,j,x,y;
3889
register MagickRealType
3890
*k,t;
3891
k=kernel->values;
3892
for( i=0, x=kernel->width-1; i<=x; i++, x--)
3893
for( j=0, y=kernel->height-1; j<y; j++, y--)
3894
{ t = k[i+j*kernel->width];
3895
k[i+j*kernel->width] = k[j+x*kernel->width];
3896
k[j+x*kernel->width] = k[x+y*kernel->width];
3897
k[x+y*kernel->width] = k[y+i*kernel->width];
3898
k[y+i*kernel->width] = t;
3899
}
3900
}
3901
/* rotate the origin - relative to center of array */
3902
{ register ssize_t x,y;
3903
x = (ssize_t) (kernel->x*2-kernel->width+1);
3904
y = (ssize_t) (kernel->y*2-kernel->height+1);
3905
kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
3906
kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
3907
}
3908
angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
3909
kernel->angle = fmod(kernel->angle+90.0, 360.0);
3910
}
3911
else
3912
perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
3913
}
3914
1850
if ( 135.0 < angle && angle <= 225.0 )
3915
1851
{
3916
/* For a 180 degree rotation - also know as a reflection
3917
* This is actually a very very common operation!
3918
* Basically all that is needed is a reversal of the kernel data!
3919
* And a reflection of the origon
3920
*/
3921
size_t
1852
/* For a 180 degree rotation - also know as a reflection */
1853
/* This is actually a very very common operation! */
1854
/* Basically all that is needed is a reversal of the kernel data! */
1855
unsigned long
3922
1856
i,j;
3923
1857
register double
3924
1858
*k,t;
3927
1861
for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
3928
1862
t=k[i], k[i]=k[j], k[j]=t;
3929
1863
3930
kernel->x = (ssize_t) kernel->width - kernel->x - 1;
3931
kernel->y = (ssize_t) kernel->height - kernel->y - 1;
3932
angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
3933
kernel->angle = fmod(kernel->angle+180.0, 360.0);
3934
}
3935
/* At this point angle should at least between -45 (315) and +45 degrees
1864
kernel->x = (long) kernel->width - kernel->x - 1;
1865
kernel->y = (long) kernel->height - kernel->y - 1;
1866
angle = fmod(angle+180.0, 360.0);
1867
}
1868
if ( 45.0 < angle && angle <= 135.0 )
1869
{ /* Do a transpose and a flop, of the image, which results in a 90
1870
* degree rotation using two mirror operations.
1871
*
1872
* WARNING: this assumes the original image was a 1 dimentional image
1873
* but currently that is the only built-ins it is applied to.
1874
*/
1875
long
1876
t;
1877
t = (long) kernel->width;
1878
kernel->width = kernel->height;
1879
kernel->height = (unsigned long) t;
1880
t = kernel->x;
1881
kernel->x = kernel->y;
1882
kernel->y = t;
1883
angle = fmod(450.0 - angle, 360.0);
1884
}
1885
/* At this point angle should be between -45 (315) and +45 degrees
3936
1886
* In the future some form of non-orthogonal angled rotates could be
3937
1887
* performed here, posibily with a linear kernel restriction.
3938
1888
*/
3939
1889
3940
return;
3941
}
3942
3943
/*
3944
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3945
% %
3946
% %
3947
% %
3948
% S c a l e G e o m e t r y K e r n e l I n f o %
3949
% %
3950
% %
3951
% %
3952
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3953
%
3954
% ScaleGeometryKernelInfo() takes a geometry argument string, typically
3955
% provided as a "-set option:convolve:scale {geometry}" user setting,
3956
% and modifies the kernel according to the parsed arguments of that setting.
3957
%
3958
% The first argument (and any normalization flags) are passed to
3959
% ScaleKernelInfo() to scale/normalize the kernel. The second argument
3960
% is then passed to UnityAddKernelInfo() to add a scled unity kernel
3961
% into the scaled/normalized kernel.
3962
%
3963
% The format of the ScaleKernelInfo method is:
3964
%
3965
% void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
3966
% const MagickStatusType normalize_flags )
3967
%
3968
% A description of each parameter follows:
3969
%
3970
% o kernel: the Morphology/Convolution kernel to modify
3971
%
3972
% o geometry:
3973
% The geometry string to parse, typically from the user provided
3974
% "-set option:convolve:scale {geometry}" setting.
3975
%
3976
*/
3977
MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
3978
const char *geometry)
3979
{
3980
GeometryFlags
3981
flags;
3982
GeometryInfo
3983
args;
3984
3985
SetGeometryInfo(&args);
3986
flags = (GeometryFlags) ParseGeometry(geometry, &args);
3987
3988
1890
#if 0
3989
/* For Debugging Geometry Input */
3990
fprintf(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
3991
flags, args.rho, args.sigma, args.xi, args.psi );
1891
Not currently in use!
1892
{ /* Do a flop, this assumes kernel is horizontally symetrical.
1893
* Each row of the kernel needs to be reversed!
1894
*/
1895
unsigned long
1896
y;
1897
register long
1898
x,r;
1899
register double
1900
*k,t;
1901
1902
for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
1903
for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
1904
t=k[x], k[x]=k[r], k[r]=t;
1905
1906
kernel->x = kernel->width - kernel->x - 1;
1907
angle = fmod(angle+180.0, 360.0);
1908
}
3992
1909
#endif
3993
3994
if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
3995
args.rho *= 0.01, args.sigma *= 0.01;
3996
3997
if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
3998
args.rho = 1.0;
3999
if ( (flags & SigmaValue) == 0 )
4000
args.sigma = 0.0;
4001
4002
/* Scale/Normalize the input kernel */
4003
ScaleKernelInfo(kernel, args.rho, flags);
4004
4005
/* Add Unity Kernel, for blending with original */
4006
if ( (flags & SigmaValue) != 0 )
4007
UnityAddKernelInfo(kernel, args.sigma);
4008
4009
1910
return;
4010
1911
}
1912
4011
1913
/*
4012
1914
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4013
1915
% %
4019
1921
% %
4020
1922
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4021
1923
%
4022
% ScaleKernelInfo() scales the given kernel list by the given amount, with or
4023
% without normalization of the sum of the kernel values (as per given flags).
1924
% ScaleKernelInfo() scales the kernel by the given amount, with or without
1925
% normalization of the sum of the kernel values.
4024
1926
%
4025
1927
% By default (no flags given) the values within the kernel is scaled
4026
% directly using given scaling factor without change.
1928
% according the given scaling factor.
4027
1929
%
4028
% If either of the two 'normalize_flags' are given the kernel will first be
4029
% normalized and then further scaled by the scaling factor value given.
1930
% If any 'normalize_flags' are given the kernel will be normalized and then
1931
% further scaled by the scaleing factor value given. A 'PercentValue' flag
1932
% will cause the given scaling factor to be divided by one hundred percent.
4030
1933
%
4031
1934
% Kernel normalization ('normalize_flags' given) is designed to ensure that
4032
1935
% any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4033
% morphology methods will fall into -1.0 to +1.0 range. Note that for
4034
% non-HDRI versions of IM this may cause images to have any negative results
4035
% clipped, unless some 'bias' is used.
1936
% morphology methods will fall into -1.0 to +1.0 range. Note however that
1937
% for non-HDRI versions of IM this may cause images to have any negative
1938
% results clipped, unless some 'clip' any negative output from 'Convolve'
1939
% with the use of some kernels.
4036
1940
%
4037
1941
% More specifically. Kernels which only contain positive values (such as a
4038
1942
% 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4039
% ensuring a 0.0 to +1.0 output range for non-HDRI images.
1943
% ensuring a 0.0 to +1.0 convolution output range for non-HDRI images.
4040
1944
%
4041
1945
% For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4042
1946
% the kernel will be scaled by the absolute of the sum of kernel values, so
4052
1956
% values seperatally to those of the negative values, so the kernel will be
4053
1957
% forced to become a zero-sum kernel better suited to such searches.
4054
1958
%
4055
% WARNING: Correct normalization of the kernel assumes that the '*_range'
1959
% WARNING: Correct normalization of the kernal assumes that the '*_range'
4056
1960
% attributes within the kernel structure have been correctly set during the
4057
1961
% kernels creation.
4058
1962
%
4059
1963
% NOTE: The values used for 'normalize_flags' have been selected specifically
4060
1964
% to match the use of geometry options, so that '!' means NormalizeValue, '^'
4061
% means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
1965
% means CorrelateNormalizeValue, and '%' means PercentValue. All other
1966
% GeometryFlags values are ignored.
4062
1967
%
4063
1968
% The format of the ScaleKernelInfo method is:
4064
1969
%
4078
1983
% specifically: NormalizeValue, CorrelateNormalizeValue,
4079
1984
% and/or PercentValue
4080
1985
%
1986
% This function is internal to this module only at this time, but can be
1987
% exported to other modules if needed.
4081
1988
*/
4082
1989
MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4083
1990
const double scaling_factor,const GeometryFlags normalize_flags)
4084
1991
{
4085
register ssize_t
1992
register long
4086
1993
i;
4087
1994
4088
1995
register double
4089
1996
pos_scale,
4090
1997
neg_scale;
4091
1998
4092
/* do the other kernels in a multi-kernel list first */
4093
if ( kernel->next != (KernelInfo *) NULL)
4094
ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4095
4096
/* Normalization of Kernel */
4097
1999
pos_scale = 1.0;
4098
2000
if ( (normalize_flags&NormalizeValue) != 0 ) {
2001
/* normalize kernel appropriately */
4099
2002
if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4100
/* non-zero-summing kernel (generally positive) */
4101
2003
pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4102
2004
else
4103
/* zero-summing kernel */
4104
pos_scale = kernel->positive_range;
2005
pos_scale = kernel->positive_range; /* special zero-summing kernel */
4105
2006
}
4106
/* Force kernel into a normalized zero-summing kernel */
2007
/* force kernel into being a normalized zero-summing kernel */
4107
2008
if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4108
2009
pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4109
2010
? kernel->positive_range : 1.0;
4116
2017
/* finialize scaling_factor for positive and negative components */
4117
2018
pos_scale = scaling_factor/pos_scale;
4118
2019
neg_scale = scaling_factor/neg_scale;
2020
if ( (normalize_flags&PercentValue) != 0 ) {
2021
pos_scale /= 100.0;
2022
neg_scale /= 100.0;
2023
}
4119
2024
4120
for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2025
for (i=0; i < (long) (kernel->width*kernel->height); i++)
4121
2026
if ( ! IsNan(kernel->values[i]) )
4122
2027
kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4123
2028
4128
2033
kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4129
2034
kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4130
2035
4131
/* swap kernel settings if user's scaling factor is negative */
2036
/* swap kernel settings if user scaling factor is negative */
4132
2037
if ( scaling_factor < MagickEpsilon ) {
4133
2038
double t;
4134
2039
t = kernel->positive_range;
4147
2052
% %
4148
2053
% %
4149
2054
% %
4150
% S h o w K e r n e l I n f o %
2055
+ S h o w K e r n e l I n f o %
4151
2056
% %
4152
2057
% %
4153
2058
% %
4164
2069
%
4165
2070
% o kernel: the Morphology/Convolution kernel
4166
2071
%
2072
% This function is internal to this module only at this time. That may change
2073
% in the future.
4167
2074
*/
4168
2075
MagickExport void ShowKernelInfo(KernelInfo *kernel)
4169
2076
{
4170
KernelInfo
4171
*k;
4172
4173
size_t
4174
c, i, u, v;
4175
4176
for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4177
4178
fprintf(stderr, "Kernel");
4179
if ( kernel->next != (KernelInfo *) NULL )
4180
fprintf(stderr, " #%lu", (unsigned long) c );
4181
fprintf(stderr, " \"%s",
4182
MagickOptionToMnemonic(MagickKernelOptions, k->type) );
4183
if ( fabs(k->angle) > MagickEpsilon )
4184
fprintf(stderr, "@%lg", k->angle);
4185
fprintf(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long) k->width,
4186
(unsigned long) k->height,(long) k->x,(long) k->y);
4187
fprintf(stderr,
4188
" with values from %.*lg to %.*lg\n",
4189
GetMagickPrecision(), k->minimum,
4190
GetMagickPrecision(), k->maximum);
4191
fprintf(stderr, "Forming a output range from %.*lg to %.*lg",
4192
GetMagickPrecision(), k->negative_range,
4193
GetMagickPrecision(), k->positive_range);
4194
if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4195
fprintf(stderr, " (Zero-Summing)\n");
4196
else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4197
fprintf(stderr, " (Normalized)\n");
4198
else
4199
fprintf(stderr, " (Sum %.*lg)\n",
4200
GetMagickPrecision(), k->positive_range+k->negative_range);
4201
for (i=v=0; v < k->height; v++) {
4202
fprintf(stderr, "%2lu:", (unsigned long) v );
4203
for (u=0; u < k->width; u++, i++)
4204
if ( IsNan(k->values[i]) )
4205
fprintf(stderr," %*s", GetMagickPrecision()+3, "nan");
4206
else
4207
fprintf(stderr," %*.*lg", GetMagickPrecision()+3,
4208
GetMagickPrecision(), k->values[i]);
4209
fprintf(stderr,"\n");
4210
}
2077
long
2078
i, u, v;
2079
2080
fprintf(stderr,
2081
"Kernel \"%s\" of size %lux%lu%+ld%+ld with values from %.*lg to %.*lg\n",
2082
MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
2083
kernel->width, kernel->height,
2084
kernel->x, kernel->y,
2085
GetMagickPrecision(), kernel->minimum,
2086
GetMagickPrecision(), kernel->maximum);
2087
fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2088
GetMagickPrecision(), kernel->negative_range,
2089
GetMagickPrecision(), kernel->positive_range,
2090
/*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2091
for (i=v=0; v < (long) kernel->height; v++) {
2092
fprintf(stderr,"%2ld:",v);
2093
for (u=0; u < (long) kernel->width; u++, i++)
2094
if ( IsNan(kernel->values[i]) )
2095
fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2096
else
2097
fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2098
GetMagickPrecision(), kernel->values[i]);
2099
fprintf(stderr,"\n");
4211
2100
}
4212
2101
}
4213
2102
4216
2105
% %
4217
2106
% %
4218
2107
% %
4219
% U n i t y A d d K e r n a l I n f o %
4220
% %
4221
% %
4222
% %
4223
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4224
%
4225
% UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4226
% to the given pre-scaled and normalized Kernel. This in effect adds that
4227
% amount of the original image into the resulting convolution kernel. This
4228
% value is usually provided by the user as a percentage value in the
4229
% 'convolve:scale' setting.
4230
%
4231
% The resulting effect is to convert the defined kernels into blended
4232
% soft-blurs, unsharp kernels or into sharpening kernels.
4233
%
4234
% The format of the UnityAdditionKernelInfo method is:
4235
%
4236
% void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4237
%
4238
% A description of each parameter follows:
4239
%
4240
% o kernel: the Morphology/Convolution kernel
4241
%
4242
% o scale:
4243
% scaling factor for the unity kernel to be added to
4244
% the given kernel.
4245
%
4246
*/
4247
MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4248
const double scale)
4249
{
4250
/* do the other kernels in a multi-kernel list first */
4251
if ( kernel->next != (KernelInfo *) NULL)
4252
UnityAddKernelInfo(kernel->next, scale);
4253
4254
/* Add the scaled unity kernel to the existing kernel */
4255
kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4256
CalcKernelMetaData(kernel); /* recalculate the meta-data */
4257
4258
return;
4259
}
4260
4261
/*
4262
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4263
% %
4264
% %
4265
% %
4266
% Z e r o K e r n e l N a n s %
2108
+ Z e r o K e r n e l N a n s %
4267
2109
% %
4268
2110
% %
4269
2111
% %
4276
2118
%
4277
2119
% The format of the ZeroKernelNans method is:
4278
2120
%
4279
% void ZeroKernelNans (KernelInfo *kernel)
2121
% voidZeroKernelNans (KernelInfo *kernel)
4280
2122
%
4281
2123
% A description of each parameter follows:
4282
2124
%
4283
2125
% o kernel: the Morphology/Convolution kernel
4284
2126
%
2127
% FUTURE: return the information in a string for API usage.
4285
2128
*/
4286
2129
MagickExport void ZeroKernelNans(KernelInfo *kernel)
4287
2130
{
4288
register size_t
2131
register long
4289
2132
i;
4290
2133
4291
/* do the other kernels in a multi-kernel list first */
4292
if ( kernel->next != (KernelInfo *) NULL)
4293
ZeroKernelNans(kernel->next);
4294
4295
for (i=0; i < (kernel->width*kernel->height); i++)
2134
for (i=0; i < (long) (kernel->width*kernel->height); i++)
4296
2135
if ( IsNan(kernel->values[i]) )
4297
2136
kernel->values[i] = 0.0;
4298
2137
Loggerhead is a web-based interface for Breezy
Version: 2.0.1