1 /*====================================================================*
2 - Copyright (C) 2001 Leptonica. All rights reserved.
3 - This software is distributed in the hope that it will be
4 - useful, but with NO WARRANTY OF ANY KIND.
5 - No author or distributor accepts responsibility to anyone for the
6 - consequences of using this software, or for whether it serves any
7 - particular purpose or works at all, unless he or she says so in
8 - writing. Everyone is granted permission to copy, modify and
9 - redistribute this source code, for commercial or non-commercial
10 - purposes, with the following restrictions: (1) the origin of this
11 - source code must not be misrepresented; (2) modified versions must
12 - be plainly marked as such; and (3) this notice may not be removed
13 - or altered from any source or modified source distribution.
14 *====================================================================*/
15
16 /*
17 * enhance.c
18 *
19 * Gamma TRC (tone reproduction curve) mapping
20 * PIX *pixGammaTRC()
21 * PIX *pixGammaTRCMasked()
22 * NUMA *numaGammaTRC()
23 *
24 * Contrast enhancement
25 * PIX *pixContrastTRC()
26 * PIX *pixContrastTRCMasked()
27 * NUMA *numaContrastTRC()
28 *
29 * Histogram equalization
30 * PIX *pixEqualizeTRC()
31 * NUMA *numaEqualizeTRC()
32 *
33 * Generic TRC mapper
34 * PIX *pixTRCMap()
35 *
36 * Unsharp-masking
37 * PIX *pixUnsharpMasking()
38 * PIX *pixUnsharpMaskingGray()
39 * PIX *pixUnsharpMaskingFast()
40 * PIX *pixUnsharpMaskingGrayFast()
41 * PIX *pixUnsharpMaskingGray1D()
42 * PIX *pixUnsharpMaskingGray2D()
43 *
44 * Hue and saturation modification
45 * PIX *pixModifyHue()
46 * PIX *pixModifySaturation()
47 * l_int32 pixMeasureSaturation()
48 *
49 * General multiplicative constant color transform
50 * PIX *pixMultConstantColor()
51 * PIX *pixMultMatrixColor()
52 *
53 * Edge by bandpass
54 * PIX *pixHalfEdgeByBandpass()
55 *
56 * Gamma correction, contrast enhancement and histogram equalization
57 * apply a simple mapping function to each pixel (or, for color
58 * images, to each sample (i.e., r,g,b) of the pixel).
59 *
60 * - Gamma correction either lightens the image or darkens
61 * it, depending on whether the gamma factor is greater
62 * or less than 1.0, respectively.
63 *
64 * - Contrast enhancement darkens the pixels that are already
65 * darker than the middle of the dynamic range (128)
66 * and lightens pixels that are lighter than 128.
67 *
68 * - Histogram equalization remaps to have the same number
69 * of image pixels at each of 256 intensity values. This is
70 * a quick and dirty method of adjusting contrast and brightness
71 * to bring out details in both light and dark regions.
72 *
73 * Unsharp masking is a more complicated enhancement.
74 * A "high frequency" image, generated by subtracting
75 * the smoothed ("low frequency") part of the image from
76 * itself, has all the energy at the edges. This "edge image"
77 * has 0 average value. A fraction of the edge image is
78 * then added to the original, enhancing the differences
79 * between pixel values at edges. Because we represent
80 * images as l_uint8 arrays, we preserve dynamic range and
81 * handle negative values by doing all the arithmetic on
82 * shifted l_uint16 arrays; the l_uint8 values are recovered
83 * at the end.
84 *
85 * Hue and saturation modification work in HSV space. Because
86 * this is too large for efficient table lookup, each pixel value
87 * is transformed to HSV, modified, and transformed back.
88 * It's not the fastest way to do this, but the method is
89 * easily understood.
90 */
91
92
93 #include <stdio.h>
94 #include <stdlib.h>
95 #include <math.h>
96 #include "allheaders.h"
97
98 /* Scales contrast enhancement factor to have a useful range
99 * between 0.0 and 1.0 */
100 static const l_float32 ENHANCE_SCALE_FACTOR = 5.;
101
102 /* Default number of pixels sampled to determine histogram */
103 static const l_int32 DEFAULT_HISTO_SAMPLES = 100000;
104
105
106 /*-------------------------------------------------------------*
107 * Gamma TRC (tone reproduction curve) mapping *
108 *-------------------------------------------------------------*/
109 /*!
110 * pixGammaTRC()
111 *
112 * Input: pixd (<optional> null or equal to pixs)
113 * pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
114 * gamma (gamma correction; must be > 0.0)
115 * minval (input value that gives 0 for output; can be < 0)
116 * maxval (input value that gives 255 for output; can be > 255)
117 * Return: pixd always
118 *
119 * Notes:
120 * (1) pixd must either be null or equal to pixs.
121 * Set pixd == pixs to get in-place operation;
122 * set pixd == null to get new image.
123 * (2) If pixs is colormapped, the colormap is transformed,
124 * either in-place or in a copy of pixs.
125 * (3) We use a gamma mapping between minval and maxval.
126 * (4) If gamma < 1.0, the image will appear darker;
127 * if gamma > 1.0, the image will appear lighter;
128 * if gamma == 1.0 and minval == 0 and maxval == 255, return a clone
129 * (5) For color images that are not colormapped, the mapping
130 * is applied to each component.
131 * (6) minval and maxval are not restricted to the interval [0, 255].
132 * If minval < 0, an input value of 0 is mapped to a
133 * nonzero output. This will turn black to gray.
134 * If maxval > 255, an input value of 255 is mapped to
135 * an output value less than 255. This will turn
136 * white (e.g., in the background) to gray.
137 * (7) Increasing minval darkens the image.
138 * (8) Decreasing maxval bleaches the image.
139 * (9) Simultaneously increasing minval and decreasing maxval
140 * will darken the image and make the colors more intense;
141 * e.g., minval = 50, maxval = 200.
142 * (10) See numaGammaTRC() for further examples of use.
143 */
144 PIX *
pixGammaTRC(PIX * pixd,PIX * pixs,l_float32 gamma,l_int32 minval,l_int32 maxval)145 pixGammaTRC(PIX *pixd,
146 PIX *pixs,
147 l_float32 gamma,
148 l_int32 minval,
149 l_int32 maxval)
150 {
151 l_int32 d;
152 NUMA *nag;
153 PIXCMAP *cmap;
154
155 PROCNAME("pixGammaTRC");
156
157 if (!pixs)
158 return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
159 if (pixd && (pixd != pixs))
160 return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
161 if (gamma <= 0.0) {
162 L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
163 gamma = 1.0;
164 }
165 if (minval >= maxval)
166 return (PIX *)ERROR_PTR("minval not < maxval", procName, pixd);
167
168 if (gamma == 1.0 && minval == 0 && maxval == 255)
169 return pixClone(pixs);
170
171 cmap = pixGetColormap(pixs);
172 d = pixGetDepth(pixs);
173 if (!cmap && d != 8 && d != 32)
174 return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
175
176 if (!pixd) /* start with a copy if not in-place */
177 pixd = pixCopy(NULL, pixs);
178
179 if (cmap) {
180 pixcmapGammaTRC(pixGetColormap(pixd), gamma, minval, maxval);
181 return pixd;
182 }
183
184 /* pixd is 8 or 32 bpp */
185 if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
186 return (PIX *)ERROR_PTR("nag not made", procName, pixd);
187 pixTRCMap(pixd, NULL, nag);
188 numaDestroy(&nag);
189
190 return pixd;
191 }
192
193
194 /*!
195 * pixGammaTRCMasked()
196 *
197 * Input: pixd (<optional> null or equal to pixs)
198 * pixs (8 or 32 bpp; not colormapped)
199 * pixm (<optional> null or 1 bpp)
200 * gamma (gamma correction; must be > 0.0)
201 * minval (input value that gives 0 for output; can be < 0)
202 * maxval (input value that gives 255 for output; can be > 255)
203 * Return: pixd always
204 *
205 * Notes:
206 * (1) Same as pixGammaTRC() except mapping is optionally over
207 * a subset of pixels described by pixm.
208 * (2) Masking does not work for colormapped images.
209 * (3) See pixGammaTRC() for details on how to use the parameters.
210 */
211 PIX *
pixGammaTRCMasked(PIX * pixd,PIX * pixs,PIX * pixm,l_float32 gamma,l_int32 minval,l_int32 maxval)212 pixGammaTRCMasked(PIX *pixd,
213 PIX *pixs,
214 PIX *pixm,
215 l_float32 gamma,
216 l_int32 minval,
217 l_int32 maxval)
218 {
219 l_int32 d;
220 NUMA *nag;
221
222 PROCNAME("pixGammaTRCMasked");
223
224 if (!pixm)
225 return pixGammaTRC(pixd, pixs, gamma, minval, maxval);
226
227 if (!pixs)
228 return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
229 if (pixGetColormap(pixs))
230 return (PIX *)ERROR_PTR("invalid: pixs has a colormap", procName, pixd);
231 if (pixd && (pixd != pixs))
232 return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
233 d = pixGetDepth(pixs);
234 if (d != 8 && d != 32)
235 return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
236 if (minval >= maxval)
237 return (PIX *)ERROR_PTR("minval not < maxval", procName, pixd);
238
239 if (gamma <= 0.0) {
240 L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
241 gamma = 1.0;
242 }
243
244 if (!pixd) /* start with a copy if not in-place */
245 pixd = pixCopy(NULL, pixs);
246
247 if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
248 return (PIX *)ERROR_PTR("nag not made", procName, pixd);
249 pixTRCMap(pixd, pixm, nag);
250 numaDestroy(&nag);
251
252 return pixd;
253 }
254
255
256 /*!
257 * numaGammaTRC()
258 *
259 * Input: gamma (gamma factor; must be > 0.0)
260 * minval (input value that gives 0 for output)
261 * maxval (input value that gives 255 for output)
262 * Return: na, or null on error
263 *
264 * Notes:
265 * (1) The map is returned as a numa; values are clipped to [0, 255].
266 * (2) To force all intensities into a range within fraction delta
267 * of white, use: minval = -256 * (1 - delta) / delta
268 * maxval = 255
269 * (3) To force all intensities into a range within fraction delta
270 * of black, use: minval = 0
271 * maxval = 256 * (1 - delta) / delta
272 */
273 NUMA *
numaGammaTRC(l_float32 gamma,l_int32 minval,l_int32 maxval)274 numaGammaTRC(l_float32 gamma,
275 l_int32 minval,
276 l_int32 maxval)
277 {
278 l_int32 i, val;
279 l_float32 x, invgamma;
280 NUMA *na;
281
282 PROCNAME("numaGammaTRC");
283
284 if (minval >= maxval)
285 return (NUMA *)ERROR_PTR("minval not < maxval", procName, NULL);
286 if (gamma <= 0.0) {
287 L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
288 gamma = 1.0;
289 }
290
291 invgamma = 1. / gamma;
292 na = numaCreate(256);
293 for (i = 0; i < minval; i++)
294 numaAddNumber(na, 0);
295 for (i = minval; i <= maxval; i++) {
296 if (i < 0) continue;
297 if (i > 255) continue;
298 x = (l_float32)(i - minval) / (l_float32)(maxval - minval);
299 val = (l_int32)(255. * powf(x, invgamma) + 0.5);
300 val = L_MAX(val, 0);
301 val = L_MIN(val, 255);
302 numaAddNumber(na, val);
303 }
304 for (i = maxval + 1; i < 256; i++)
305 numaAddNumber(na, 255);
306
307 return na;
308 }
309
310
311 /*-------------------------------------------------------------*
312 * Contrast enhancement *
313 *-------------------------------------------------------------*/
314 /*!
315 * pixContrastTRC()
316 *
317 * Input: pixd (<optional> null or equal to pixs)
318 * pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
319 * factor (0.0 is no enhancement)
320 * Return: pixd always
321 *
322 * Notes:
323 * (1) pixd must either be null or equal to pixs.
324 * Set pixd == pixs to get in-place operation;
325 * set pixd == null to get new image.
326 * (2) If pixs is colormapped, the colormap is transformed,
327 * either in-place or in a copy of pixs.
328 * (3) Contrast is enhanced by mapping each color component
329 * using an atan function with maximum slope at 127.
330 * Pixels below 127 are lowered in intensity and pixels
331 * above 127 are increased.
332 * (4) The useful range for the contrast factor is scaled to
333 * be in (0.0 to 1.0), but larger values can also be used.
334 * 0.0 corresponds to no enhancement; return a clone.
335 * (5) For color images that are not colormapped, the mapping
336 * is applied to each component.
337 */
338 PIX *
pixContrastTRC(PIX * pixd,PIX * pixs,l_float32 factor)339 pixContrastTRC(PIX *pixd,
340 PIX *pixs,
341 l_float32 factor)
342 {
343 l_int32 d;
344 NUMA *nac;
345 PIXCMAP *cmap;
346
347 PROCNAME("pixContrastTRC");
348
349 if (!pixs)
350 return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
351 if (pixd && (pixd != pixs))
352 return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
353 if (factor < 0.0) {
354 L_WARNING("factor must be >= 0.0; using 0.0", procName);
355 factor = 0.0;
356 }
357 if (factor == 0.0)
358 return pixClone(pixs);
359
360 cmap = pixGetColormap(pixs);
361 d = pixGetDepth(pixs);
362 if (!cmap && d != 8 && d != 32)
363 return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
364
365 if (!pixd) /* start with a copy if not in-place */
366 pixd = pixCopy(NULL, pixs);
367
368 if (cmap) {
369 pixcmapContrastTRC(pixGetColormap(pixd), factor);
370 return pixd;
371 }
372
373 /* pixd is 8 or 32 bpp */
374 if ((nac = numaContrastTRC(factor)) == NULL)
375 return (PIX *)ERROR_PTR("nac not made", procName, pixd);
376 pixTRCMap(pixd, NULL, nac);
377 numaDestroy(&nac);
378
379 return pixd;
380 }
381
382
383 /*!
384 * pixContrastTRCMasked()
385 *
386 * Input: pixd (<optional> null or equal to pixs)
387 * pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
388 * pixm (<optional> null or 1 bpp)
389 * factor (0.0 is no enhancement)
390 * Return: pixd always
391 *
392 * Notes:
393 * (1) Same as pixContrastTRC() except mapping is optionally over
394 * a subset of pixels described by pixm.
395 * (2) Masking does not work for colormapped images.
396 * (3) See pixContrastTRC() for details on how to use the parameters.
397 */
398 PIX *
pixContrastTRCMasked(PIX * pixd,PIX * pixs,PIX * pixm,l_float32 factor)399 pixContrastTRCMasked(PIX *pixd,
400 PIX *pixs,
401 PIX *pixm,
402 l_float32 factor)
403 {
404 l_int32 d;
405 NUMA *nac;
406
407 PROCNAME("pixContrastTRCMasked");
408
409 if (!pixm)
410 return pixContrastTRC(pixd, pixs, factor);
411
412 if (!pixs)
413 return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
414 if (pixGetColormap(pixs))
415 return (PIX *)ERROR_PTR("invalid: pixs has a colormap", procName, pixd);
416 if (pixd && (pixd != pixs))
417 return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
418 d = pixGetDepth(pixs);
419 if (d != 8 && d != 32)
420 return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
421
422 if (factor < 0.0) {
423 L_WARNING("factor must be >= 0.0; using 0.0", procName);
424 factor = 0.0;
425 }
426 if (factor == 0.0)
427 return pixClone(pixs);
428
429 if (!pixd) /* start with a copy if not in-place */
430 pixd = pixCopy(NULL, pixs);
431
432 if ((nac = numaContrastTRC(factor)) == NULL)
433 return (PIX *)ERROR_PTR("nac not made", procName, pixd);
434 pixTRCMap(pixd, pixm, nac);
435 numaDestroy(&nac);
436
437 return pixd;
438 }
439
440
441 /*!
442 * numaContrastTRC()
443 *
444 * Input: factor (generally between 0.0 (no enhancement)
445 * and 1.0, but can be larger than 1.0)
446 * Return: na, or null on error
447 *
448 * Notes:
449 * (1) The mapping is monotonic increasing, where 0 is mapped
450 * to 0 and 255 is mapped to 255.
451 * (2) As 'factor' is increased from 0.0 (where the mapping is linear),
452 * the map gets closer to its limit as a step function that
453 * jumps from 0 to 255 at the center (input value = 127).
454 */
455 NUMA *
numaContrastTRC(l_float32 factor)456 numaContrastTRC(l_float32 factor)
457 {
458 l_int32 i, val;
459 l_float64 x, ymax, ymin, dely, scale;
460 NUMA *na;
461
462 PROCNAME("numaContrastTRC");
463
464 if (factor < 0.0) {
465 L_WARNING("factor must be >= 0.0; using 0.0; no enhancement", procName);
466 factor = 0.0;
467 }
468 if (factor == 0.0)
469 return numaMakeSequence(0, 1, 256); /* linear map */
470
471 scale = ENHANCE_SCALE_FACTOR;
472 ymax = atan((l_float64)(1.0 * factor * scale));
473 ymin = atan((l_float64)(-127. * factor * scale / 128.));
474 dely = ymax - ymin;
475 na = numaCreate(256);
476 for (i = 0; i < 256; i++) {
477 x = (l_float64)i;
478 val = (l_int32)((255. / dely) *
479 (-ymin + atan((l_float64)(factor * scale * (x - 127.) / 128.))) +
480 0.5);
481 numaAddNumber(na, val);
482 }
483
484 return na;
485 }
486
487
488 /*-------------------------------------------------------------*
489 * Histogram equalization *
490 *-------------------------------------------------------------*/
491 /*!
492 * pixEqualizeTRC()
493 *
494 * Input: pixd (<optional> null or equal to pixs)
495 * pixs (8 bpp, or colormapped)
496 * fract (fraction of equalization movement of pixel values)
497 * factor (subsampling factor; integer >= 1)
498 * Return: pixd, or null on error
499 *
500 * Notes:
501 * (1) pixd must either be null or equal to pixs.
502 * Set pixd == pixs to get in-place operation;
503 * set pixd == null to get new image.
504 * (2) In histogram equalization, a tone reproduction curve
505 * mapping is used to make the number of pixels at each
506 * intensity equal.
507 * (3) If fract == 0.0, no equalization is performed; return a copy.
508 * If fract == 1.0, equalization is complete.
509 * (4) Set the subsampling factor > 1 to reduce the amount of computation.
510 * (5) If the pix is colormapped, the colormap is removed and
511 * it is converted to grayscale.
512 * (6) Note that even if there is a colormap, we can get an
513 * in-place operation because the intermediate image pixt
514 * is copied back to pixs (which for in-place is the same
515 * as pixd).
516 */
517 PIX *
pixEqualizeTRC(PIX * pixd,PIX * pixs,l_float32 fract,l_int32 factor)518 pixEqualizeTRC(PIX *pixd,
519 PIX *pixs,
520 l_float32 fract,
521 l_int32 factor)
522 {
523 NUMA *na;
524 PIX *pixt;
525 PIXCMAP *cmap;
526
527 PROCNAME("pixEqualizeTRC");
528
529 if (!pixs)
530 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
531 if (pixd && (pixd != pixs))
532 return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
533 cmap = pixGetColormap(pixs);
534 if ((pixGetDepth(pixs) != 8) && !cmap)
535 return (PIX *)ERROR_PTR("pixs not 8 bpp or cmapped", procName, NULL);
536 if (fract < 0.0 || fract > 1.0)
537 return (PIX *)ERROR_PTR("fract not in [0.0 ... 1.0]", procName, NULL);
538 if (factor < 1)
539 return (PIX *)ERROR_PTR("sampling factor < 1", procName, NULL);
540
541 if (fract == 0.0)
542 return pixCopy(pixd, pixs);
543
544 /* If there is a colormap, remove it. */
545 if (cmap)
546 pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
547 else
548 pixt = pixClone(pixs);
549
550 /* Make a copy if necessary */
551 pixd = pixCopy(pixd, pixt);
552 pixDestroy(&pixt);
553
554 if ((na = numaEqualizeTRC(pixd, fract, factor)) == NULL)
555 return (PIX *)ERROR_PTR("na not made", procName, pixd);
556 pixTRCMap(pixd, NULL, na);
557 numaDestroy(&na);
558
559 return pixd;
560 }
561
562
563 /*!
564 * numaEqualizeTRC()
565 *
566 * Input: pix (8 bpp, no colormap)
567 * fract (fraction of equalization movement of pixel values)
568 * factor (subsampling factor; integer >= 1)
569 * Return: nad, or null on error
570 *
571 * Notes:
572 * (1) If fract == 0.0, no equalization will be performed.
573 * If fract == 1.0, equalization is complete.
574 * (2) Set the subsampling factor > 1 to reduce the amount of computation.
575 * (3) The map is returned as a numa with 256 values, specifying
576 * the equalized value (array value) for every input value
577 * (the array index).
578 */
579 NUMA *
numaEqualizeTRC(PIX * pix,l_float32 fract,l_int32 factor)580 numaEqualizeTRC(PIX *pix,
581 l_float32 fract,
582 l_int32 factor)
583 {
584 l_int32 iin, iout, itarg;
585 l_float32 val, sum;
586 NUMA *nah, *nasum, *nad;
587
588 PROCNAME("numaEqualizeTRC");
589
590 if (!pix)
591 return (NUMA *)ERROR_PTR("pix not defined", procName, NULL);
592 if (pixGetDepth(pix) != 8)
593 return (NUMA *)ERROR_PTR("pix not 8 bpp", procName, NULL);
594 if (fract < 0.0 || fract > 1.0)
595 return (NUMA *)ERROR_PTR("fract not in [0.0 ... 1.0]", procName, NULL);
596 if (factor < 1)
597 return (NUMA *)ERROR_PTR("sampling factor < 1", procName, NULL);
598
599 if (fract == 0.0)
600 L_WARNING("fract = 0.0; no equalization requested", procName);
601
602 if ((nah = pixGetGrayHistogram(pix, factor)) == NULL)
603 return (NUMA *)ERROR_PTR("histogram not made", procName, NULL);
604 numaGetSum(nah, &sum);
605 nasum = numaGetPartialSums(nah);
606
607 nad = numaCreate(256);
608 for (iin = 0; iin < 256; iin++) {
609 numaGetFValue(nasum, iin, &val);
610 itarg = (l_int32)(255. * val / sum + 0.5);
611 iout = iin + (l_int32)(fract * (itarg - iin));
612 iout = L_MIN(iout, 255); /* to be safe */
613 numaAddNumber(nad, iout);
614 }
615
616 numaDestroy(&nah);
617 numaDestroy(&nasum);
618 return nad;
619 }
620
621
622 /*-------------------------------------------------------------*
623 * Generic TRC mapping *
624 *-------------------------------------------------------------*/
625 /*!
626 * pixTRCMap()
627 *
628 * Input: pixs (8 grayscale or 32 bpp rgb; not colormapped)
629 * pixm (<optional> 1 bpp mask)
630 * na (mapping array)
631 * Return: pixd, or null on error
632 *
633 * Notes:
634 * (1) This operation is in-place on pixs.
635 * (2) For 32 bpp, this applies the same map to each of the r,g,b
636 * components.
637 * (3) The mapping array is of size 256, and it maps the input
638 * index into values in the range [0, 255].
639 * (4) If defined, the optional 1 bpp mask pixm has its origin
640 * aligned with pixs, and the map function is applied only
641 * to pixels in pixs under the fg of pixm.
642 */
643 l_int32
pixTRCMap(PIX * pixs,PIX * pixm,NUMA * na)644 pixTRCMap(PIX *pixs,
645 PIX *pixm,
646 NUMA *na)
647 {
648 l_int32 w, h, d, wm, hm, wpl, wplm, i, j, sval8, dval8;
649 l_int32 *tab;
650 l_uint32 sval32, dval32;
651 l_uint32 *data, *datam, *line, *linem;
652
653 PROCNAME("pixTRCMap");
654
655 if (!pixs)
656 return ERROR_INT("pixs not defined", procName, 1);
657 if (pixGetColormap(pixs))
658 return ERROR_INT("pixs is colormapped", procName, 1);
659 if (!na)
660 return ERROR_INT("na not defined", procName, 1);
661 if (numaGetCount(na) != 256)
662 return ERROR_INT("na not of size 256", procName, 1);
663 pixGetDimensions(pixs, &w, &h, &d);
664 if (d != 8 && d != 32)
665 return ERROR_INT("pixs not 8 or 32 bpp", procName, 1);
666 if (pixm) {
667 if (pixGetDepth(pixm) != 1)
668 return ERROR_INT("pixm not 1 bpp", procName, 1);
669 }
670
671 tab = numaGetIArray(na); /* get the array for efficiency */
672 wpl = pixGetWpl(pixs);
673 data = pixGetData(pixs);
674 if (!pixm) {
675 if (d == 8) {
676 for (i = 0; i < h; i++) {
677 line = data + i * wpl;
678 for (j = 0; j < w; j++) {
679 sval8 = GET_DATA_BYTE(line, j);
680 dval8 = tab[sval8];
681 SET_DATA_BYTE(line, j, dval8);
682 }
683 }
684 }
685 else { /* d == 32 */
686 for (i = 0; i < h; i++) {
687 line = data + i * wpl;
688 for (j = 0; j < w; j++) {
689 sval32 = *(line + j);
690 dval32 =
691 tab[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
692 tab[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
693 tab[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
694 *(line + j) = dval32;
695 }
696 }
697 }
698 }
699 else {
700 datam = pixGetData(pixm);
701 wplm = pixGetWpl(pixm);
702 pixGetDimensions(pixm, &wm, &hm, NULL);
703 if (d == 8) {
704 for (i = 0; i < h; i++) {
705 if (i >= hm)
706 break;
707 line = data + i * wpl;
708 linem = datam + i * wplm;
709 for (j = 0; j < w; j++) {
710 if (j >= wm)
711 break;
712 if (GET_DATA_BIT(linem, j) == 0)
713 continue;
714 sval8 = GET_DATA_BYTE(line, j);
715 dval8 = tab[sval8];
716 SET_DATA_BYTE(line, j, dval8);
717 }
718 }
719 }
720 else { /* d == 32 */
721 for (i = 0; i < h; i++) {
722 if (i >= hm)
723 break;
724 line = data + i * wpl;
725 linem = datam + i * wplm;
726 for (j = 0; j < w; j++) {
727 if (j >= wm)
728 break;
729 if (GET_DATA_BIT(linem, j) == 0)
730 continue;
731 sval32 = *(line + j);
732 dval32 =
733 tab[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
734 tab[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
735 tab[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
736 *(line + j) = dval32;
737 }
738 }
739 }
740 }
741
742 FREE(tab);
743 return 0;
744 }
745
746
747
748 /*-----------------------------------------------------------------------*
749 * Unsharp masking *
750 *-----------------------------------------------------------------------*/
751 /*!
752 * pixUnsharpMasking()
753 *
754 * Input: pixs (all depths except 1 bpp; with or without colormaps)
755 * halfwidth ("half-width" of smoothing filter)
756 * fract (fraction of edge added back into image)
757 * Return: pixd, or null on error
758 *
759 * Notes:
760 * (1) We use symmetric smoothing filters of odd dimension,
761 * typically use sizes of 3, 5, 7, etc. The @halfwidth parameter
762 * for these is (size - 1)/2; i.e., 1, 2, 3, etc.
763 * (2) The fract parameter is typically taken in the
764 * range: 0.2 < fract < 0.7
765 */
766 PIX *
pixUnsharpMasking(PIX * pixs,l_int32 halfwidth,l_float32 fract)767 pixUnsharpMasking(PIX *pixs,
768 l_int32 halfwidth,
769 l_float32 fract)
770 {
771 l_int32 d;
772 PIX *pixt, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;
773
774 PROCNAME("pixUnsharpMasking");
775
776 if (!pixs || (pixGetDepth(pixs) == 1))
777 return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", procName, NULL);
778 if (fract <= 0.0 || halfwidth <= 0) {
779 L_WARNING("no sharpening requested; clone returned", procName);
780 return pixClone(pixs);
781 }
782
783 if (halfwidth == 1 || halfwidth == 2)
784 return pixUnsharpMaskingFast(pixs, halfwidth, fract, L_BOTH_DIRECTIONS);
785
786 /* Remove colormap; clone if possible; result is either 8 or 32 bpp */
787 if ((pixt = pixConvertTo8Or32(pixs, 0, 1)) == NULL)
788 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
789
790 /* Sharpen */
791 d = pixGetDepth(pixt);
792 if (d == 8)
793 pixd = pixUnsharpMaskingGray(pixt, halfwidth, fract);
794 else { /* d == 32 */
795 pixr = pixGetRGBComponent(pixs, COLOR_RED);
796 pixrs = pixUnsharpMaskingGray(pixr, halfwidth, fract);
797 pixDestroy(&pixr);
798 pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
799 pixgs = pixUnsharpMaskingGray(pixg, halfwidth, fract);
800 pixDestroy(&pixg);
801 pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
802 pixbs = pixUnsharpMaskingGray(pixb, halfwidth, fract);
803 pixDestroy(&pixb);
804 pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
805 pixDestroy(&pixrs);
806 pixDestroy(&pixgs);
807 pixDestroy(&pixbs);
808 }
809
810 pixDestroy(&pixt);
811 return pixd;
812 }
813
814
815 /*!
816 * pixUnsharpMaskingGray()
817 *
818 * Input: pixs (8 bpp; no colormap)
819 * halfwidth ("half-width" of smoothing filter)
820 * fract (fraction of edge added back into image)
821 * Return: pixd, or null on error
822 *
823 * Notes:
824 * (1) We use symmetric smoothing filters of odd dimension,
825 * typically use sizes of 3, 5, 7, etc. The @halfwidth parameter
826 * for these is (size - 1)/2; i.e., 1, 2, 3, etc.
827 * (2) The fract parameter is typically taken in the range:
828 * 0.2 < fract < 0.7
829 */
830 PIX *
pixUnsharpMaskingGray(PIX * pixs,l_int32 halfwidth,l_float32 fract)831 pixUnsharpMaskingGray(PIX *pixs,
832 l_int32 halfwidth,
833 l_float32 fract)
834 {
835 l_int32 w, h, d;
836 PIX *pixc, *pixd;
837 PIXACC *pixacc;
838
839 PROCNAME("pixUnsharpMaskingGray");
840
841 if (!pixs)
842 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
843 pixGetDimensions(pixs, &w, &h, &d);
844 if (d != 8 || pixGetColormap(pixs) != NULL)
845 return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
846 if (fract <= 0.0 || halfwidth <= 0) {
847 L_WARNING("no sharpening requested; clone returned", procName);
848 return pixClone(pixs);
849 }
850 if (halfwidth == 1 || halfwidth == 2)
851 return pixUnsharpMaskingGrayFast(pixs, halfwidth, fract,
852 L_BOTH_DIRECTIONS);
853
854 if ((pixc = pixBlockconvGray(pixs, NULL, halfwidth, halfwidth)) == NULL)
855 return (PIX *)ERROR_PTR("pixc not made", procName, NULL);
856
857 /* Steps:
858 * (1) edge image is pixs - pixc (this is highpass part)
859 * (2) multiply edge image by fract
860 * (3) add fraction of edge to pixs
861 *
862 * To show how this is done with both interfaces to arithmetic
863 * on integer Pix, here is the implementation in the lower-level
864 * function calls:
865 * pixt = pixInitAccumulate(w, h, 0x10000000)) == NULL)
866 * pixAccumulate(pixt, pixs, L_ARITH_ADD);
867 * pixAccumulate(pixt, pixc, L_ARITH_SUBTRACT);
868 * pixMultConstAccumulate(pixt, fract, 0x10000000);
869 * pixAccumulate(pixt, pixs, L_ARITH_ADD);
870 * pixd = pixFinalAccumulate(pixt, 0x10000000, 8)) == NULL)
871 * pixDestroy(&pixt);
872 *
873 * The code below does the same thing using the Pixacc accumulator,
874 * hiding the details of the offset that is needed for subtraction.
875 */
876 pixacc = pixaccCreate(w, h, 1);
877 pixaccAdd(pixacc, pixs);
878 pixaccSubtract(pixacc, pixc);
879 pixaccMultConst(pixacc, fract);
880 pixaccAdd(pixacc, pixs);
881 pixd = pixaccFinal(pixacc, 8);
882 pixaccDestroy(&pixacc);
883
884 pixDestroy(&pixc);
885 return pixd;
886 }
887
888
889 /*!
890 * pixUnsharpMaskingFast()
891 *
892 * Input: pixs (all depths except 1 bpp; with or without colormaps)
893 * halfwidth ("half-width" of smoothing filter; 1 and 2 only)
894 * fract (fraction of high frequency added to image)
895 * direction (L_HORIZ, L_VERT, L_BOTH_DIRECTIONS)
896 * Return: pixd, or null on error
897 *
898 * Notes:
899 * (1) The fast version uses separable 1-D filters directly on
900 * the input image. The halfwidth is either 1 (full width = 3)
901 * or 2 (full width = 5).
902 * (2) The fract parameter is typically taken in the
903 * range: 0.2 < fract < 0.7
904 * (3) To skip horizontal sharpening, use @fracth = 0.0; ditto for @fractv
905 * (4) For one dimensional filtering (as an example):
906 * For @halfwidth = 1, the low-pass filter is
907 * L: 1/3 1/3 1/3
908 * and the high-pass filter is
909 * H = I - L: -1/3 2/3 -1/3
910 * For @halfwidth = 2, the low-pass filter is
911 * L: 1/5 1/5 1/5 1/5 1/5
912 * and the high-pass filter is
913 * H = I - L: -1/5 -1/5 4/5 -1/5 -1/5
914 * The new sharpened pixel value is found by adding some fraction
915 * of the high-pass filter value (which sums to 0) to the
916 * initial pixel value:
917 * N = I + fract * H
918 * (5) For 2D, the sharpening filter is not separable, because the
919 * vertical filter depends on the horizontal location relative
920 * to the filter origin, and v.v. So we either do the full
921 * 2D filter (for @halfwidth == 1) or do the low-pass
922 * convolution separably and then compose with the original pix.
923 */
924 PIX *
pixUnsharpMaskingFast(PIX * pixs,l_int32 halfwidth,l_float32 fract,l_int32 direction)925 pixUnsharpMaskingFast(PIX *pixs,
926 l_int32 halfwidth,
927 l_float32 fract,
928 l_int32 direction)
929 {
930 l_int32 d;
931 PIX *pixt, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;
932
933 PROCNAME("pixUnsharpMaskingFast");
934
935 if (!pixs || (pixGetDepth(pixs) == 1))
936 return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", procName, NULL);
937 if (fract <= 0.0 || halfwidth <= 0) {
938 L_WARNING("no sharpening requested; clone returned", procName);
939 return pixClone(pixs);
940 }
941 if (halfwidth != 1 && halfwidth != 2)
942 return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
943 if (direction != L_HORIZ && direction != L_VERT &&
944 direction != L_BOTH_DIRECTIONS)
945 return (PIX *)ERROR_PTR("invalid direction", procName, NULL);
946
947 /* Remove colormap; clone if possible; result is either 8 or 32 bpp */
948 if ((pixt = pixConvertTo8Or32(pixs, 0, 1)) == NULL)
949 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
950
951 /* Sharpen */
952 d = pixGetDepth(pixt);
953 if (d == 8)
954 pixd = pixUnsharpMaskingGrayFast(pixt, halfwidth, fract, direction);
955 else { /* d == 32 */
956 pixr = pixGetRGBComponent(pixs, COLOR_RED);
957 pixrs = pixUnsharpMaskingGrayFast(pixr, halfwidth, fract, direction);
958 pixDestroy(&pixr);
959 pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
960 pixgs = pixUnsharpMaskingGrayFast(pixg, halfwidth, fract, direction);
961 pixDestroy(&pixg);
962 pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
963 pixbs = pixUnsharpMaskingGrayFast(pixb, halfwidth, fract, direction);
964 pixDestroy(&pixb);
965 pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
966 pixDestroy(&pixrs);
967 pixDestroy(&pixgs);
968 pixDestroy(&pixbs);
969 }
970
971 pixDestroy(&pixt);
972 return pixd;
973 }
974
975
976
977 /*!
978 * pixUnsharpMaskingGrayFast()
979 *
980 * Input: pixs (8 bpp; no colormap)
981 * halfwidth ("half-width" of smoothing filter: 1 or 2)
982 * fract (fraction of high frequency added to image)
983 * direction (L_HORIZ, L_VERT, L_BOTH_DIRECTIONS)
984 * Return: pixd, or null on error
985 *
986 * Notes:
987 * (1) For usage and explanation of the algorithm, see notes
988 * in pixUnsharpMaskingFast().
989 */
990 PIX *
pixUnsharpMaskingGrayFast(PIX * pixs,l_int32 halfwidth,l_float32 fract,l_int32 direction)991 pixUnsharpMaskingGrayFast(PIX *pixs,
992 l_int32 halfwidth,
993 l_float32 fract,
994 l_int32 direction)
995 {
996 PIX *pixd;
997
998 PROCNAME("pixUnsharpMaskingGrayFast");
999
1000 if (!pixs)
1001 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1002 if (pixGetDepth(pixs) != 8 || pixGetColormap(pixs) != NULL)
1003 return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
1004 if (fract <= 0.0 || halfwidth <= 0) {
1005 L_WARNING("no sharpening requested; clone returned", procName);
1006 return pixClone(pixs);
1007 }
1008 if (halfwidth != 1 && halfwidth != 2)
1009 return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
1010 if (direction != L_HORIZ && direction != L_VERT &&
1011 direction != L_BOTH_DIRECTIONS)
1012 return (PIX *)ERROR_PTR("invalid direction", procName, NULL);
1013
1014 if (direction != L_BOTH_DIRECTIONS)
1015 pixd = pixUnsharpMaskingGray1D(pixs, halfwidth, fract, direction);
1016 else /* 2D sharpening */
1017 pixd = pixUnsharpMaskingGray2D(pixs, halfwidth, fract);
1018
1019 return pixd;
1020 }
1021
1022
1023 /*!
1024 * pixUnsharpMaskingGray1D()
1025 *
1026 * Input: pixs (8 bpp; no colormap)
1027 * halfwidth ("half-width" of smoothing filter: 1 or 2)
1028 * fract (fraction of high frequency added to image)
1029 * direction (of filtering; use L_HORIZ or L_VERT)
1030 * Return: pixd, or null on error
1031 *
1032 * Notes:
1033 * (1) For usage and explanation of the algorithm, see notes
1034 * in pixUnsharpMaskingFast().
1035 */
1036 PIX *
pixUnsharpMaskingGray1D(PIX * pixs,l_int32 halfwidth,l_float32 fract,l_int32 direction)1037 pixUnsharpMaskingGray1D(PIX *pixs,
1038 l_int32 halfwidth,
1039 l_float32 fract,
1040 l_int32 direction)
1041 {
1042 l_int32 w, h, d, wpls, wpld, i, j, ival;
1043 l_uint32 *datas, *datad;
1044 l_uint32 *lines, *lines0, *lines1, *lines2, *lines3, *lines4, *lined;
1045 l_float32 val, a[5];
1046 PIX *pixd;
1047
1048 PROCNAME("pixUnsharpMaskingGray1D");
1049
1050 if (!pixs)
1051 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1052 pixGetDimensions(pixs, &w, &h, &d);
1053 if (d != 8 || pixGetColormap(pixs) != NULL)
1054 return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
1055 if (fract <= 0.0 || halfwidth <= 0) {
1056 L_WARNING("no sharpening requested; clone returned", procName);
1057 return pixClone(pixs);
1058 }
1059 if (halfwidth != 1 && halfwidth != 2)
1060 return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
1061
1062 /* Initialize pixd with pixels from pixs that will not be
1063 * set when computing the sharpened values. */
1064 pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
1065 halfwidth, halfwidth);
1066 datas = pixGetData(pixs);
1067 datad = pixGetData(pixd);
1068 wpls = pixGetWpl(pixs);
1069 wpld = pixGetWpl(pixd);
1070
1071 if (halfwidth == 1) {
1072 a[0] = -fract / 3.0;
1073 a[1] = 1.0 + fract * 2.0 / 3.0;
1074 a[2] = a[0];
1075 }
1076 else { /* halfwidth == 2 */
1077 a[0] = -fract / 5.0;
1078 a[1] = a[0];
1079 a[2] = 1.0 + fract * 4.0 / 5.0;
1080 a[3] = a[0];
1081 a[4] = a[0];
1082 }
1083
1084 if (direction == L_HORIZ) {
1085 for (i = 0; i < h; i++) {
1086 lines = datas + i * wpls;
1087 lined = datad + i * wpld;
1088 if (halfwidth == 1) {
1089 for (j = 1; j < w - 1; j++) {
1090 val = a[0] * GET_DATA_BYTE(lines, j - 1) +
1091 a[1] * GET_DATA_BYTE(lines, j) +
1092 a[2] * GET_DATA_BYTE(lines, j + 1);
1093 ival = (l_int32)val;
1094 ival = L_MAX(0, ival);
1095 ival = L_MIN(255, ival);
1096 SET_DATA_BYTE(lined, j, ival);
1097 }
1098 }
1099 else { /* halfwidth == 2 */
1100 for (j = 2; j < w - 2; j++) {
1101 val = a[0] * GET_DATA_BYTE(lines, j - 2) +
1102 a[1] * GET_DATA_BYTE(lines, j - 1) +
1103 a[2] * GET_DATA_BYTE(lines, j) +
1104 a[3] * GET_DATA_BYTE(lines, j + 1) +
1105 a[4] * GET_DATA_BYTE(lines, j + 2);
1106 ival = (l_int32)val;
1107 ival = L_MAX(0, ival);
1108 ival = L_MIN(255, ival);
1109 SET_DATA_BYTE(lined, j, ival);
1110 }
1111 }
1112 }
1113 }
1114 else { /* direction == L_VERT */
1115 if (halfwidth == 1) {
1116 for (i = 1; i < h - 1; i++) {
1117 lines0 = datas + (i - 1) * wpls;
1118 lines1 = datas + i * wpls;
1119 lines2 = datas + (i + 1) * wpls;
1120 lined = datad + i * wpld;
1121 for (j = 0; j < w; j++) {
1122 val = a[0] * GET_DATA_BYTE(lines0, j) +
1123 a[1] * GET_DATA_BYTE(lines1, j) +
1124 a[2] * GET_DATA_BYTE(lines2, j);
1125 ival = (l_int32)val;
1126 ival = L_MAX(0, ival);
1127 ival = L_MIN(255, ival);
1128 SET_DATA_BYTE(lined, j, ival);
1129 }
1130 }
1131 }
1132 else { /* halfwidth == 2 */
1133 for (i = 2; i < h - 2; i++) {
1134 lines0 = datas + (i - 2) * wpls;
1135 lines1 = datas + (i - 1) * wpls;
1136 lines2 = datas + i * wpls;
1137 lines3 = datas + (i + 1) * wpls;
1138 lines4 = datas + (i + 2) * wpls;
1139 lined = datad + i * wpld;
1140 for (j = 0; j < w; j++) {
1141 val = a[0] * GET_DATA_BYTE(lines0, j) +
1142 a[1] * GET_DATA_BYTE(lines1, j) +
1143 a[2] * GET_DATA_BYTE(lines2, j) +
1144 a[3] * GET_DATA_BYTE(lines3, j) +
1145 a[4] * GET_DATA_BYTE(lines4, j);
1146 ival = (l_int32)val;
1147 ival = L_MAX(0, ival);
1148 ival = L_MIN(255, ival);
1149 SET_DATA_BYTE(lined, j, ival);
1150 }
1151 }
1152 }
1153 }
1154
1155 return pixd;
1156 }
1157
1158
1159 /*!
1160 * pixUnsharpMaskingGray2D()
1161 *
1162 * Input: pixs (8 bpp; no colormap)
1163 * halfwidth ("half-width" of smoothing filter: 1 or 2)
1164 * fract (fraction of high frequency added to image)
1165 * Return: pixd, or null on error
1166 *
1167 * Notes:
1168 * (1) For halfwidth == 1, we implement the full sharpening filter
1169 * directly. For halfwidth == 2, we implement the the lowpass
1170 * filter separably and then compute the sharpening result locally.
1171 */
1172 PIX *
pixUnsharpMaskingGray2D(PIX * pixs,l_int32 halfwidth,l_float32 fract)1173 pixUnsharpMaskingGray2D(PIX *pixs,
1174 l_int32 halfwidth,
1175 l_float32 fract)
1176 {
1177 l_int32 w, h, d, wpls, wpld, wplf, i, j, ival, sval;
1178 l_uint32 *datas, *datad, *lines, *lines0, *lines1, *lines2, *lined;
1179 l_float32 val, a[9];
1180 l_float32 *dataf, *linef, *linef0, *linef1, *linef2, *linef3, *linef4;
1181 PIX *pixd;
1182 FPIX *fpix;
1183
1184 PROCNAME("pixUnsharpMaskingGray2D");
1185
1186 if (!pixs)
1187 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1188 pixGetDimensions(pixs, &w, &h, &d);
1189 if (d != 8 || pixGetColormap(pixs) != NULL)
1190 return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
1191 if (fract <= 0.0 || halfwidth <= 0) {
1192 L_WARNING("no sharpening requested; clone returned", procName);
1193 return pixClone(pixs);
1194 }
1195 if (halfwidth != 1 && halfwidth != 2)
1196 return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
1197
1198 pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
1199 halfwidth, halfwidth);
1200 datad = pixGetData(pixd);
1201 wpld = pixGetWpl(pixd);
1202 datas = pixGetData(pixs);
1203 wpls = pixGetWpl(pixs);
1204
1205 if (halfwidth == 1) {
1206 for (i = 0; i < 9; i++)
1207 a[i] = -fract / 9.0;
1208 a[4] = 1.0 + fract * 8.0 / 9.0;
1209 for (i = 1; i < h - 1; i++) {
1210 lines0 = datas + (i - 1) * wpls;
1211 lines1 = datas + i * wpls;
1212 lines2 = datas + (i + 1) * wpls;
1213 lined = datad + i * wpld;
1214 for (j = 1; j < w - 1; j++) {
1215 val = a[0] * GET_DATA_BYTE(lines0, j - 1) +
1216 a[1] * GET_DATA_BYTE(lines0, j) +
1217 a[2] * GET_DATA_BYTE(lines0, j + 1) +
1218 a[3] * GET_DATA_BYTE(lines1, j - 1) +
1219 a[4] * GET_DATA_BYTE(lines1, j) +
1220 a[5] * GET_DATA_BYTE(lines1, j + 1) +
1221 a[6] * GET_DATA_BYTE(lines2, j - 1) +
1222 a[7] * GET_DATA_BYTE(lines2, j) +
1223 a[8] * GET_DATA_BYTE(lines2, j + 1);
1224 ival = (l_int32)(val + 0.5);
1225 ival = L_MAX(0, ival);
1226 ival = L_MIN(255, ival);
1227 SET_DATA_BYTE(lined, j, ival);
1228 }
1229 }
1230
1231 return pixd;
1232 }
1233
1234 /* For halfwidth == 2, do the low pass separably. Store
1235 * the result of horizontal smoothing in an intermediate fpix. */
1236 fpix = fpixCreate(w, h);
1237 dataf = fpixGetData(fpix);
1238 wplf = fpixGetWpl(fpix);
1239 for (i = 2; i < h - 2; i++) {
1240 lines = datas + i * wpls;
1241 linef = dataf + i * wplf;
1242 for (j = 2; j < w - 2; j++) {
1243 val = GET_DATA_BYTE(lines, j - 2) +
1244 GET_DATA_BYTE(lines, j - 1) +
1245 GET_DATA_BYTE(lines, j) +
1246 GET_DATA_BYTE(lines, j + 1) +
1247 GET_DATA_BYTE(lines, j + 2);
1248 linef[j] = val;
1249 }
1250 }
1251
1252 /* Do vertical smoothing to finish the low-pass filter.
1253 * At each pixel, if L is the lowpass value, I is the
1254 * src pixel value and f is the fraction of highpass to
1255 * be added to I, then the highpass filter value is
1256 * H = I - L
1257 * and the new sharpened value is
1258 * N = I + f * H.
1259 */
1260 for (i = 2; i < h - 2; i++) {
1261 linef0 = dataf + (i - 2) * wplf;
1262 linef1 = dataf + (i - 1) * wplf;
1263 linef2 = dataf + i * wplf;
1264 linef3 = dataf + (i + 1) * wplf;
1265 linef4 = dataf + (i + 2) * wplf;
1266 lined = datad + i * wpld;
1267 lines = datas + i * wpls;
1268 for (j = 2; j < w - 2; j++) {
1269 val = 0.04 * (linef0[j] + linef1[j] + linef2[j] +
1270 linef3[j] + linef4[j]); /* L: lowpass filter value */
1271 sval = GET_DATA_BYTE(lines, j); /* I: source pixel */
1272 ival = (l_int32)(sval + fract * (sval - val) + 0.5);
1273 ival = L_MAX(0, ival);
1274 ival = L_MIN(255, ival);
1275 SET_DATA_BYTE(lined, j, ival);
1276 }
1277 }
1278
1279 fpixDestroy(&fpix);
1280 return pixd;
1281 }
1282
1283
1284
1285 /*-----------------------------------------------------------------------*
1286 * Hue and saturation modification *
1287 *-----------------------------------------------------------------------*/
1288 /*!
1289 * pixModifyHue()
1290 *
1291 * Input: pixd (<optional> can be null, existing or equal to pixs)
1292 * pixs (32 bpp rgb)
1293 * fract (between -1.0 and 1.0)
1294 * Return: pixd, or null on error
1295 *
1296 * Notes:
1297 * (1) Use fract > 0.0 to increase hue value; < 0.0 to decrease it.
1298 * 1.0 (or -1.0) represents a 360 degree rotation; i.e., no change.
1299 */
1300 PIX *
pixModifyHue(PIX * pixd,PIX * pixs,l_float32 fract)1301 pixModifyHue(PIX *pixd,
1302 PIX *pixs,
1303 l_float32 fract)
1304 {
1305 l_int32 w, h, d, i, j, wpl, delhue;
1306 l_int32 rval, gval, bval, hval, sval, vval;
1307 l_uint32 *data, *line;
1308
1309 PROCNAME("pixModifyHue");
1310
1311 if (!pixs)
1312 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1313 pixGetDimensions(pixs, &w, &h, &d);
1314 if (d != 32)
1315 return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
1316 if (L_ABS(fract) > 1.0)
1317 return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", procName, NULL);
1318
1319 pixd = pixCopy(pixd, pixs);
1320
1321 delhue = (l_int32)(240 * fract);
1322 if (delhue == 0 || delhue == 240 || delhue == -240) {
1323 L_WARNING("no change requested in hue", procName);
1324 return pixd;
1325 }
1326 if (delhue < 0)
1327 delhue += 240;
1328
1329 data = pixGetData(pixd);
1330 wpl = pixGetWpl(pixd);
1331 for (i = 0; i < h; i++) {
1332 line = data + i * wpl;
1333 for (j = 0; j < w; j++) {
1334 extractRGBValues(line[j], &rval, &gval, &bval);
1335 convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
1336 hval = (hval + delhue) % 240;
1337 convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
1338 composeRGBPixel(rval, gval, bval, line + j);
1339 }
1340 }
1341
1342 return pixd;
1343 }
1344
1345
1346 /*!
1347 * pixModifySaturation()
1348 *
1349 * Input: pixd (<optional> can be null, existing or equal to pixs)
1350 * pixs (32 bpp rgb)
1351 * fract (between -1.0 and 1.0)
1352 * Return: pixd, or null on error
1353 *
1354 * Notes:
1355 * (1) If fract > 0.0, it gives the fraction that the pixel
1356 * saturation is moved from its initial value toward 255.
1357 * If fract < 0.0, it gives the fraction that the pixel
1358 * saturation is moved from its initial value toward 0.
1359 * The limiting values for fract = -1.0 (1.0) thus set the
1360 * saturation to 0 (255).
1361 */
1362 PIX *
pixModifySaturation(PIX * pixd,PIX * pixs,l_float32 fract)1363 pixModifySaturation(PIX *pixd,
1364 PIX *pixs,
1365 l_float32 fract)
1366 {
1367 l_int32 w, h, d, i, j, wpl;
1368 l_int32 rval, gval, bval, hval, sval, vval;
1369 l_uint32 *data, *line;
1370
1371 PROCNAME("pixModifySaturation");
1372
1373 if (!pixs)
1374 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1375 pixGetDimensions(pixs, &w, &h, &d);
1376 if (d != 32)
1377 return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
1378 if (L_ABS(fract) > 1.0)
1379 return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", procName, NULL);
1380
1381 pixd = pixCopy(pixd, pixs);
1382 if (fract == 0.0) {
1383 L_WARNING("no change requested in saturation", procName);
1384 return pixd;
1385 }
1386
1387 data = pixGetData(pixd);
1388 wpl = pixGetWpl(pixd);
1389 for (i = 0; i < h; i++) {
1390 line = data + i * wpl;
1391 for (j = 0; j < w; j++) {
1392 extractRGBValues(line[j], &rval, &gval, &bval);
1393 convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
1394 if (fract < 0.0)
1395 sval = (l_int32)(sval * (1.0 + fract));
1396 else
1397 sval = (l_int32)(sval + fract * (255 - sval));
1398 convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
1399 composeRGBPixel(rval, gval, bval, line + j);
1400 }
1401 }
1402
1403 return pixd;
1404 }
1405
1406
1407 /*!
1408 * pixMeasureSaturation()
1409 *
1410 * Input: pixs (32 bpp rgb)
1411 * factor (subsampling factor; integer >= 1)
1412 * &sat (<return> average saturation)
1413 * Return: pixd, or null on error
1414 */
1415 l_int32
pixMeasureSaturation(PIX * pixs,l_int32 factor,l_float32 * psat)1416 pixMeasureSaturation(PIX *pixs,
1417 l_int32 factor,
1418 l_float32 *psat)
1419 {
1420 l_int32 w, h, d, i, j, wpl, sum, count;
1421 l_int32 rval, gval, bval, hval, sval, vval;
1422 l_uint32 *data, *line;
1423
1424 PROCNAME("pixMeasureSaturation");
1425
1426 if (!psat)
1427 return ERROR_INT("pixs not defined", procName, 1);
1428 *psat = 0.0;
1429 if (!pixs)
1430 return ERROR_INT("pixs not defined", procName, 1);
1431 pixGetDimensions(pixs, &w, &h, &d);
1432 if (d != 32)
1433 return ERROR_INT("pixs not 32 bpp", procName, 1);
1434 if (factor < 1)
1435 return ERROR_INT("subsampling factor < 1", procName, 1);
1436
1437 data = pixGetData(pixs);
1438 wpl = pixGetWpl(pixs);
1439 for (i = 0, sum = 0, count = 0; i < h; i += factor) {
1440 line = data + i * wpl;
1441 for (j = 0; j < w; j += factor) {
1442 extractRGBValues(line[j], &rval, &gval, &bval);
1443 convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
1444 sum += sval;
1445 count++;
1446 }
1447 }
1448
1449 if (count > 0)
1450 *psat = (l_float32)sum / (l_float32)count;
1451 return 0;
1452 }
1453
1454
1455 /*-----------------------------------------------------------------------*
1456 * General multiplicative constant color transform *
1457 *-----------------------------------------------------------------------*/
1458 /*
1459 * pixMultConstantColor()
1460 *
1461 * Input: pixs (colormapped or rgb)
1462 * rfact, gfact, bfact (multiplicative factors on each component)
1463 * Return: pixd (colormapped or rgb, with colors scaled), or null on error
1464 *
1465 * Notes:
1466 * (1) rfact, gfact and bfact can only have non-negative values.
1467 * They can be greater than 1.0. All transformed component
1468 * values are clipped to the interval [0, 255].
1469 * (2) For multiplication with a general 3x3 matrix of constants,
1470 * use pixMultMatrixColor().
1471 */
1472 PIX *
pixMultConstantColor(PIX * pixs,l_float32 rfact,l_float32 gfact,l_float32 bfact)1473 pixMultConstantColor(PIX *pixs,
1474 l_float32 rfact,
1475 l_float32 gfact,
1476 l_float32 bfact)
1477 {
1478 l_int32 i, j, w, h, d, wpls, wpld;
1479 l_int32 ncolors, rval, gval, bval, nrval, ngval, nbval;
1480 l_uint32 nval;
1481 l_uint32 *datas, *datad, *lines, *lined;
1482 PIX *pixd;
1483 PIXCMAP *cmap;
1484
1485 PROCNAME("pixMultConstantColor");
1486
1487 if (!pixs)
1488 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1489 pixGetDimensions(pixs, &w, &h, &d);
1490 cmap = pixGetColormap(pixs);
1491 if (!cmap && d != 32)
1492 return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
1493 rfact = L_MAX(0.0, rfact);
1494 gfact = L_MAX(0.0, gfact);
1495 bfact = L_MAX(0.0, bfact);
1496
1497 if (cmap) {
1498 if ((pixd = pixCopy(NULL, pixs)) == NULL)
1499 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
1500 cmap = pixGetColormap(pixd);
1501 ncolors = pixcmapGetCount(cmap);
1502 for (i = 0; i < ncolors; i++) {
1503 pixcmapGetColor(cmap, i, &rval, &gval, &bval);
1504 nrval = (l_int32)(rfact * rval);
1505 ngval = (l_int32)(gfact * gval);
1506 nbval = (l_int32)(bfact * bval);
1507 nrval = L_MIN(255, nrval);
1508 ngval = L_MIN(255, ngval);
1509 nbval = L_MIN(255, nbval);
1510 pixcmapResetColor(cmap, i, nrval, ngval, nbval);
1511 }
1512 return pixd;
1513 }
1514
1515 if ((pixd = pixCreateTemplateNoInit(pixs)) == NULL)
1516 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
1517 datas = pixGetData(pixs);
1518 datad = pixGetData(pixd);
1519 wpls = pixGetWpl(pixs);
1520 wpld = pixGetWpl(pixd);
1521 for (i = 0; i < h; i++) {
1522 lines = datas + i * wpls;
1523 lined = datad + i * wpld;
1524 for (j = 0; j < w; j++) {
1525 extractRGBValues(lines[j], &rval, &gval, &bval);
1526 nrval = (l_int32)(rfact * rval);
1527 ngval = (l_int32)(gfact * gval);
1528 nbval = (l_int32)(bfact * bval);
1529 nrval = L_MIN(255, nrval);
1530 ngval = L_MIN(255, ngval);
1531 nbval = L_MIN(255, nbval);
1532 composeRGBPixel(nrval, ngval, nbval, &nval);
1533 *(lined + j) = nval;
1534 }
1535 }
1536
1537 return pixd;
1538 }
1539
1540
1541 /*
1542 * pixMultMatrixColor()
1543 *
1544 * Input: pixs (colormapped or rgb)
1545 * kernel (3x3 matrix of floats)
1546 * Return: pixd (colormapped or rgb), or null on error
1547 *
1548 * Notes:
1549 * (1) The kernel is a data structure used mostly for floating point
1550 * convolution. Here it is a 3x3 matrix of floats that are used
1551 * to transform the pixel values by matrix multiplication:
1552 * nrval = a[0,0] * rval + a[0,1] * gval + a[0,2] * bval
1553 * ngval = a[1,0] * rval + a[1,1] * gval + a[1,2] * bval
1554 * nbval = a[2,0] * rval + a[2,1] * gval + a[2,2] * bval
1555 * (2) The matrix can be generated in several ways.
1556 * See kernel.c for details. Here are two of them:
1557 * (a) kel = kernelCreate(3, 3);
1558 * kernelSetElement(kel, 0, 0, val00);
1559 * kernelSetElement(kel, 0, 1, val01);
1560 * ...
1561 * (b) from a static string; e.g.,:
1562 * const char *kdata = " 0.6 0.3 -0.2 "
1563 * " 0.1 1.2 0.4 "
1564 * " -0.4 0.2 0.9 ";
1565 * kel = kernelCreateFromString(3, 3, 0, 0, kdata);
1566 * (3) For the special case where the matrix is diagonal, it is easier
1567 * to use pixMultConstantColor().
1568 * (4) Matrix entries can have positive and negative values, and can
1569 * be larger than 1.0. All transformed component values
1570 * are clipped to [0, 255].
1571 */
1572 PIX *
pixMultMatrixColor(PIX * pixs,L_KERNEL * kel)1573 pixMultMatrixColor(PIX *pixs,
1574 L_KERNEL *kel)
1575 {
1576 l_int32 i, j, index, kw, kh, w, h, d, wpls, wpld;
1577 l_int32 ncolors, rval, gval, bval, nrval, ngval, nbval;
1578 l_uint32 nval;
1579 l_uint32 *datas, *datad, *lines, *lined;
1580 l_float32 v[9]; /* use linear array for convenience */
1581 PIX *pixd;
1582 PIXCMAP *cmap;
1583
1584 PROCNAME("pixMultMatrixColor");
1585
1586 if (!pixs)
1587 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1588 if (!kel)
1589 return (PIX *)ERROR_PTR("kel not defined", procName, NULL);
1590 kernelGetParameters(kel, &kw, &kh, NULL, NULL);
1591 if (kw != 3 || kh != 3)
1592 return (PIX *)ERROR_PTR("matrix not 3x3", procName, NULL);
1593 pixGetDimensions(pixs, &w, &h, &d);
1594 cmap = pixGetColormap(pixs);
1595 if (!cmap && d != 32)
1596 return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
1597
1598 for (i = 0, index = 0; i < 3; i++)
1599 for (j = 0; j < 3; j++, index++)
1600 kernelGetElement(kel, i, j, v + index);
1601
1602 if (cmap) {
1603 if ((pixd = pixCopy(NULL, pixs)) == NULL)
1604 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
1605 cmap = pixGetColormap(pixd);
1606 ncolors = pixcmapGetCount(cmap);
1607 for (i = 0; i < ncolors; i++) {
1608 pixcmapGetColor(cmap, i, &rval, &gval, &bval);
1609 nrval = (l_int32)(v[0] * rval + v[1] * gval + v[2] * bval);
1610 ngval = (l_int32)(v[3] * rval + v[4] * gval + v[5] * bval);
1611 nbval = (l_int32)(v[6] * rval + v[7] * gval + v[8] * bval);
1612 nrval = L_MAX(0, L_MIN(255, nrval));
1613 ngval = L_MAX(0, L_MIN(255, ngval));
1614 nbval = L_MAX(0, L_MIN(255, nbval));
1615 pixcmapResetColor(cmap, i, nrval, ngval, nbval);
1616 }
1617 return pixd;
1618 }
1619
1620 if ((pixd = pixCreateTemplateNoInit(pixs)) == NULL)
1621 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
1622 datas = pixGetData(pixs);
1623 datad = pixGetData(pixd);
1624 wpls = pixGetWpl(pixs);
1625 wpld = pixGetWpl(pixd);
1626 for (i = 0; i < h; i++) {
1627 lines = datas + i * wpls;
1628 lined = datad + i * wpld;
1629 for (j = 0; j < w; j++) {
1630 extractRGBValues(lines[j], &rval, &gval, &bval);
1631 nrval = (l_int32)(v[0] * rval + v[1] * gval + v[2] * bval);
1632 ngval = (l_int32)(v[3] * rval + v[4] * gval + v[5] * bval);
1633 nbval = (l_int32)(v[6] * rval + v[7] * gval + v[8] * bval);
1634 nrval = L_MAX(0, L_MIN(255, nrval));
1635 ngval = L_MAX(0, L_MIN(255, ngval));
1636 nbval = L_MAX(0, L_MIN(255, nbval));
1637 composeRGBPixel(nrval, ngval, nbval, &nval);
1638 *(lined + j) = nval;
1639 }
1640 }
1641
1642 return pixd;
1643 }
1644
1645
1646 /*-------------------------------------------------------------*
1647 * Half-edge by bandpass *
1648 *-------------------------------------------------------------*/
1649 /*!
1650 * pixHalfEdgeByBandpass()
1651 *
1652 * Input: pixs (8 bpp gray or 32 bpp rgb)
1653 * sm1h, sm1v ("half-widths" of smoothing filter sm1)
1654 * sm2h, sm2v ("half-widths" of smoothing filter sm2)
1655 * (require sm2 != sm1)
1656 * Return: pixd, or null on error
1657 *
1658 * Notes:
1659 * (1) We use symmetric smoothing filters of odd dimension,
1660 * typically use 3, 5, 7, etc. The smoothing parameters
1661 * for these are 1, 2, 3, etc. The filter size is related
1662 * to the smoothing parameter by
1663 * size = 2 * smoothing + 1
1664 * (2) Because we take the difference of two lowpass filters,
1665 * this is actually a bandpass filter.
1666 * (3) We allow both filters to be anisotropic.
1667 * (4) Consider either the h or v component of the 2 filters.
1668 * Depending on whether sm1 > sm2 or sm2 > sm1, we get
1669 * different halves of the smoothed gradients (or "edges").
1670 * This difference of smoothed signals looks more like
1671 * a second derivative of a transition, which we rectify
1672 * by not allowing the signal to go below zero. If sm1 < sm2,
1673 * the sm2 transition is broader, so the difference between
1674 * sm1 and sm2 signals is positive on the upper half of
1675 * the transition. Likewise, if sm1 > sm2, the sm1 - sm2
1676 * signal difference is positive on the lower half of
1677 * the transition.
1678 */
1679 PIX *
pixHalfEdgeByBandpass(PIX * pixs,l_int32 sm1h,l_int32 sm1v,l_int32 sm2h,l_int32 sm2v)1680 pixHalfEdgeByBandpass(PIX *pixs,
1681 l_int32 sm1h,
1682 l_int32 sm1v,
1683 l_int32 sm2h,
1684 l_int32 sm2v)
1685 {
1686 l_int32 d;
1687 PIX *pixg, *pixacc, *pixc1, *pixc2;
1688
1689 PROCNAME("pixHalfEdgeByBandpass");
1690
1691 if (!pixs)
1692 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1693 if (sm1h == sm2h && sm1v == sm2v)
1694 return (PIX *)ERROR_PTR("sm2 = sm1", procName, NULL);
1695 d = pixGetDepth(pixs);
1696 if (d != 8 && d != 32)
1697 return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
1698 if (d == 32)
1699 pixg = pixConvertRGBToLuminance(pixs);
1700 else /* d == 8 */
1701 pixg = pixClone(pixs);
1702
1703 /* Make a convolution accumulator and use it twice */
1704 if ((pixacc = pixBlockconvAccum(pixg)) == NULL)
1705 return (PIX *)ERROR_PTR("pixacc not made", procName, NULL);
1706 if ((pixc1 = pixBlockconvGray(pixg, pixacc, sm1h, sm1v)) == NULL)
1707 return (PIX *)ERROR_PTR("pixc1 not made", procName, NULL);
1708 if ((pixc2 = pixBlockconvGray(pixg, pixacc, sm2h, sm2v)) == NULL)
1709 return (PIX *)ERROR_PTR("pixc2 not made", procName, NULL);
1710 pixDestroy(&pixacc);
1711
1712 /* Compute the half-edge using pixc1 - pixc2. */
1713 pixSubtractGray(pixc1, pixc1, pixc2);
1714
1715 pixDestroy(&pixg);
1716 pixDestroy(&pixc2);
1717 return pixc1;
1718 }
1719
1720
1721