Coverage Report

Created: 2024-07-27 06:27

/src/leptonica/src/enhance.c
Line
Count
Source (jump to first uncovered line)
1
/*====================================================================*
2
 -  Copyright (C) 2001 Leptonica.  All rights reserved.
3
 -
4
 -  Redistribution and use in source and binary forms, with or without
5
 -  modification, are permitted provided that the following conditions
6
 -  are met:
7
 -  1. Redistributions of source code must retain the above copyright
8
 -     notice, this list of conditions and the following disclaimer.
9
 -  2. Redistributions in binary form must reproduce the above
10
 -     copyright notice, this list of conditions and the following
11
 -     disclaimer in the documentation and/or other materials
12
 -     provided with the distribution.
13
 -
14
 -  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
15
 -  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
16
 -  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
17
 -  A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL ANY
18
 -  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
19
 -  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
20
 -  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
21
 -  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
22
 -  OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
23
 -  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24
 -  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25
 *====================================================================*/
26
27
/*!
28
 * \file enhance.c
29
 * <pre>
30
 *
31
 *      Gamma TRC (tone reproduction curve) mapping
32
 *           PIX     *pixGammaTRC()
33
 *           PIX     *pixGammaTRCMasked()
34
 *           PIX     *pixGammaTRCWithAlpha()
35
 *           NUMA    *numaGammaTRC()
36
 *
37
 *      Contrast enhancement
38
 *           PIX     *pixContrastTRC()
39
 *           PIX     *pixContrastTRCMasked()
40
 *           NUMA    *numaContrastTRC()
41
 *
42
 *      Histogram equalization
43
 *           PIX     *pixEqualizeTRC()
44
 *           NUMA    *numaEqualizeTRC()
45
 *
46
 *      Generic TRC mapper
47
 *           l_int32  pixTRCMap()
48
 *           l_int32  pixTRCMapGeneral()
49
 *
50
 *      Unsharp-masking
51
 *           PIX     *pixUnsharpMasking()
52
 *           PIX     *pixUnsharpMaskingGray()
53
 *           PIX     *pixUnsharpMaskingFast()
54
 *           PIX     *pixUnsharpMaskingGrayFast()
55
 *           PIX     *pixUnsharpMaskingGray1D()
56
 *           PIX     *pixUnsharpMaskingGray2D()
57
 *
58
 *      Hue and saturation modification
59
 *           PIX     *pixModifyHue()
60
 *           PIX     *pixModifySaturation()
61
 *           l_int32  pixMeasureSaturation()
62
 *           PIX     *pixModifyBrightness()
63
 *
64
 *      Color shifting
65
 *           PIX     *pixMosaicColorShiftRGB()
66
 *           PIX     *pixColorShiftRGB()
67
 *
68
 *      Darken gray (unsaturated) pixels
69
 *           PIX     *pixDarkenGray()
70
 *
71
 *      General multiplicative constant color transform
72
 *           PIX     *pixMultConstantColor()
73
 *           PIX     *pixMultMatrixColor()
74
 *
75
 *      Edge by bandpass
76
 *           PIX     *pixHalfEdgeByBandpass()
77
 *
78
 *      Gamma correction, contrast enhancement and histogram equalization
79
 *      apply a simple mapping function to each pixel (or, for color
80
 *      images, to each sample (i.e., r,g,b) of the pixel).
81
 *
82
 *       ~ Gamma correction either lightens the image or darkens
83
 *         it, depending on whether the gamma factor is greater
84
 *         or less than 1.0, respectively.
85
 *
86
 *       ~ Contrast enhancement darkens the pixels that are already
87
 *         darker than the middle of the dynamic range (128)
88
 *         and lightens pixels that are lighter than 128.
89
 *
90
 *       ~ Histogram equalization remaps to have the same number
91
 *         of image pixels at each of 256 intensity values.  This is
92
 *         a quick and dirty method of adjusting contrast and brightness
93
 *         to bring out details in both light and dark regions.
94
 *
95
 *      Unsharp masking is a more complicated enhancement.
96
 *      A "high frequency" image, generated by subtracting
97
 *      the smoothed ("low frequency") part of the image from
98
 *      itself, has all the energy at the edges.  This "edge image"
99
 *      has 0 average value.  A fraction of the edge image is
100
 *      then added to the original, enhancing the differences
101
 *      between pixel values at edges.  Because we represent
102
 *      images as l_uint8 arrays, we preserve dynamic range and
103
 *      handle negative values by doing all the arithmetic on
104
 *      shifted l_uint16 arrays; the l_uint8 values are recovered
105
 *      at the end.
106
 *
107
 *      Hue and saturation modification work in HSV space.  Because
108
 *      this is too large for efficient table lookup, each pixel value
109
 *      is transformed to HSV, modified, and transformed back.
110
 *      It's not the fastest way to do this, but the method is
111
 *      easily understood.
112
 *
113
 *      Unsharp masking is never in-place, and returns a clone if no
114
 *      operation is to be performed.
115
 * </pre>
116
 */
117
118
#ifdef HAVE_CONFIG_H
119
#include <config_auto.h>
120
#endif  /* HAVE_CONFIG_H */
121
122
#include <math.h>
123
#include "allheaders.h"
124
125
    /* Scales contrast enhancement factor to have a useful range
126
     * between 0.0 and 1.0 */
127
static const l_float32  EnhanceScaleFactor = 5.0;
128
129
/*-------------------------------------------------------------*
130
 *         Gamma TRC (tone reproduction curve) mapping         *
131
 *-------------------------------------------------------------*/
132
/*!
133
 * \brief   pixGammaTRC()
134
 *
135
 * \param[in]    pixd     [optional] null or equal to pixs
136
 * \param[in]    pixs     8 or 32 bpp; or 2, 4 or 8 bpp with colormap
137
 * \param[in]    gamma    gamma correction; must be > 0.0
138
 * \param[in]    minval   input value that gives 0 for output; can be < 0
139
 * \param[in]    maxval   input value that gives 255 for output; can be > 255
140
 * \return  pixd always
141
 *
142
 * <pre>
143
 * Notes:
144
 *      (1) pixd must either be null or equal to pixs.
145
 *          For in-place operation, set pixd == pixs:
146
 *             pixGammaTRC(pixs, pixs, ...);
147
 *          To get a new image, set pixd == null:
148
 *             pixd = pixGammaTRC(NULL, pixs, ...);
149
 *      (2) If pixs is colormapped, the colormap is transformed,
150
 *          either in-place or in a copy of pixs.
151
 *      (3) We use a gamma mapping between minval and maxval.
152
 *      (4) If gamma < 1.0, the image will appear darker;
153
 *          if gamma > 1.0, the image will appear lighter;
154
 *      (5) If gamma = 1.0 and minval = 0 and maxval = 255, no
155
 *          enhancement is performed; return a copy unless in-place,
156
 *          in which case this is a no-op.
157
 *      (6) For color images that are not colormapped, the mapping
158
 *          is applied to each component.
159
 *      (7) minval and maxval are not restricted to the interval [0, 255].
160
 *          If minval < 0, an input value of 0 is mapped to a
161
 *          nonzero output.  This will turn black to gray.
162
 *          If maxval > 255, an input value of 255 is mapped to
163
 *          an output value less than 255.  This will turn
164
 *          white (e.g., in the background) to gray.
165
 *      (8) Increasing minval darkens the image.
166
 *      (9) Decreasing maxval bleaches the image.
167
 *      (10) Simultaneously increasing minval and decreasing maxval
168
 *           will darken the image and make the colors more intense;
169
 *           e.g., minval = 50, maxval = 200.
170
 *      (11) See numaGammaTRC() for further examples of use.
171
 *      (12) Use pixTRCMapGeneral() if applying different mappings
172
 *           to each channel in an RGB image.
173
 * </pre>
174
 */
175
PIX *
176
pixGammaTRC(PIX       *pixd,
177
            PIX       *pixs,
178
            l_float32  gamma,
179
            l_int32    minval,
180
            l_int32    maxval)
181
0
{
182
0
l_int32   d;
183
0
NUMA     *nag;
184
0
PIXCMAP  *cmap;
185
186
0
    if (!pixs)
187
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
188
0
    if (pixd && (pixd != pixs))
189
0
        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
190
0
    if (gamma <= 0.0) {
191
0
        L_WARNING("gamma must be > 0.0; setting to 1.0\n", __func__);
192
0
        gamma = 1.0;
193
0
    }
194
0
    if (minval >= maxval)
195
0
        return (PIX *)ERROR_PTR("minval not < maxval", __func__, pixd);
196
0
    cmap = pixGetColormap(pixs);
197
0
    d = pixGetDepth(pixs);
198
0
    if (!cmap && d != 8 && d != 32)
199
0
        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, pixd);
200
201
0
    if (gamma == 1.0 && minval == 0 && maxval == 255)  /* no-op */
202
0
        return pixCopy(pixd, pixs);
203
204
0
    if (!pixd)  /* start with a copy if not in-place */
205
0
        pixd = pixCopy(NULL, pixs);
206
207
0
    if (cmap) {
208
0
        pixcmapGammaTRC(pixGetColormap(pixd), gamma, minval, maxval);
209
0
        return pixd;
210
0
    }
211
212
        /* pixd is 8 or 32 bpp */
213
0
    if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
214
0
        return (PIX *)ERROR_PTR("nag not made", __func__, pixd);
215
0
    pixTRCMap(pixd, NULL, nag);
216
0
    numaDestroy(&nag);
217
218
0
    return pixd;
219
0
}
220
221
222
/*!
223
 * \brief   pixGammaTRCMasked()
224
 *
225
 * \param[in]    pixd      [optional] null or equal to pixs
226
 * \param[in]    pixs      8 or 32 bpp; not colormapped
227
 * \param[in]    pixm      [optional] null or 1 bpp
228
 * \param[in]    gamma     gamma correction; must be > 0.0
229
 * \param[in]    minval    input value that gives 0 for output; can be < 0
230
 * \param[in]    maxval    input value that gives 255 for output; can be > 255
231
 * \return  pixd always
232
 *
233
 * <pre>
234
 * Notes:
235
 *      (1) Same as pixGammaTRC() except mapping is optionally over
236
 *          a subset of pixels described by pixm.
237
 *      (2) Masking does not work for colormapped images.
238
 *      (3) See pixGammaTRC() for details on how to use the parameters.
239
 * </pre>
240
 */
241
PIX *
242
pixGammaTRCMasked(PIX       *pixd,
243
                  PIX       *pixs,
244
                  PIX       *pixm,
245
                  l_float32  gamma,
246
                  l_int32    minval,
247
                  l_int32    maxval)
248
0
{
249
0
l_int32  d;
250
0
NUMA    *nag;
251
252
0
    if (!pixm)
253
0
        return pixGammaTRC(pixd, pixs, gamma, minval, maxval);
254
255
0
    if (!pixs)
256
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
257
0
    if (pixGetColormap(pixs))
258
0
        return (PIX *)ERROR_PTR("invalid: pixs has a colormap", __func__, pixd);
259
0
    if (pixd && (pixd != pixs))
260
0
        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
261
0
    d = pixGetDepth(pixs);
262
0
    if (d != 8 && d != 32)
263
0
        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, pixd);
264
0
    if (minval >= maxval)
265
0
        return (PIX *)ERROR_PTR("minval not < maxval", __func__, pixd);
266
0
    if (gamma <= 0.0) {
267
0
        L_WARNING("gamma must be > 0.0; setting to 1.0\n", __func__);
268
0
        gamma = 1.0;
269
0
    }
270
271
0
    if (gamma == 1.0 && minval == 0 && maxval == 255)
272
0
        return pixCopy(pixd, pixs);
273
274
0
    if (!pixd)  /* start with a copy if not in-place */
275
0
        pixd = pixCopy(NULL, pixs);
276
277
0
    if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
278
0
        return (PIX *)ERROR_PTR("nag not made", __func__, pixd);
279
0
    pixTRCMap(pixd, pixm, nag);
280
0
    numaDestroy(&nag);
281
282
0
    return pixd;
283
0
}
284
285
286
/*!
287
 * \brief   pixGammaTRCWithAlpha()
288
 *
289
 * \param[in]    pixd     [optional] null or equal to pixs
290
 * \param[in]    pixs     32 bpp
291
 * \param[in]    gamma    gamma correction; must be > 0.0
292
 * \param[in]    minval   input value that gives 0 for output; can be < 0
293
 * \param[in]    maxval   input value that gives 255 for output; can be > 255
294
 * \return  pixd always
295
 *
296
 * <pre>
297
 * Notes:
298
 *      (1) See usage notes in pixGammaTRC().
299
 *      (2) This version saves the alpha channel.  It is only valid
300
 *          for 32 bpp (no colormap), and is a bit slower.
301
 * </pre>
302
 */
303
PIX *
304
pixGammaTRCWithAlpha(PIX       *pixd,
305
                     PIX       *pixs,
306
                     l_float32  gamma,
307
                     l_int32    minval,
308
                     l_int32    maxval)
309
0
{
310
0
NUMA  *nag;
311
0
PIX   *pixalpha;
312
313
0
    if (!pixs || pixGetDepth(pixs) != 32)
314
0
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", __func__, pixd);
315
0
    if (pixd && (pixd != pixs))
316
0
        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
317
0
    if (gamma <= 0.0) {
318
0
        L_WARNING("gamma must be > 0.0; setting to 1.0\n", __func__);
319
0
        gamma = 1.0;
320
0
    }
321
0
    if (minval >= maxval)
322
0
        return (PIX *)ERROR_PTR("minval not < maxval", __func__, pixd);
323
324
0
    if (gamma == 1.0 && minval == 0 && maxval == 255)
325
0
        return pixCopy(pixd, pixs);
326
0
    if (!pixd)  /* start with a copy if not in-place */
327
0
        pixd = pixCopy(NULL, pixs);
328
329
0
    pixalpha = pixGetRGBComponent(pixs, L_ALPHA_CHANNEL);  /* save */
330
0
    if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
331
0
        return (PIX *)ERROR_PTR("nag not made", __func__, pixd);
332
0
    pixTRCMap(pixd, NULL, nag);
333
0
    pixSetRGBComponent(pixd, pixalpha, L_ALPHA_CHANNEL);  /* restore */
334
0
    pixSetSpp(pixd, 4);
335
336
0
    numaDestroy(&nag);
337
0
    pixDestroy(&pixalpha);
338
0
    return pixd;
339
0
}
340
341
342
/*!
343
 * \brief   numaGammaTRC()
344
 *
345
 * \param[in]    gamma   gamma factor; must be > 0.0
346
 * \param[in]    minval  input value that gives 0 for output
347
 * \param[in]    maxval  input value that gives 255 for output
348
 * \return  na, or NULL on error
349
 *
350
 * <pre>
351
 * Notes:
352
 *      (1) The map is returned as a numa; values are clipped to [0, 255].
353
 *      (2) For a linear mapping, set gamma = 1.0.
354
 *      (3) To force all intensities into a range within fraction delta
355
 *          of white, use: minval = -256 * (1 - delta) / delta
356
 *                         maxval = 255
357
 *      (4) To force all intensities into a range within fraction delta
358
 *          of black, use: minval = 0
359
 *                         maxval = 256 * (1 - delta) / delta
360
 * </pre>
361
 */
362
NUMA *
363
numaGammaTRC(l_float32  gamma,
364
             l_int32    minval,
365
             l_int32    maxval)
366
0
{
367
0
l_int32    i, val;
368
0
l_float32  x, invgamma;
369
0
NUMA      *na;
370
371
0
    if (minval >= maxval)
372
0
        return (NUMA *)ERROR_PTR("minval not < maxval", __func__, NULL);
373
0
    if (gamma <= 0.0) {
374
0
        L_WARNING("gamma must be > 0.0; setting to 1.0\n", __func__);
375
0
        gamma = 1.0;
376
0
    }
377
378
0
    invgamma = 1. / gamma;
379
0
    na = numaCreate(256);
380
0
    for (i = 0; i < minval; i++)
381
0
        numaAddNumber(na, 0);
382
0
    for (i = minval; i <= maxval; i++) {
383
0
        if (i < 0) continue;
384
0
        if (i > 255) continue;
385
0
        x = (l_float32)(i - minval) / (l_float32)(maxval - minval);
386
0
        val = (l_int32)(255. * powf(x, invgamma) + 0.5);
387
0
        val = L_MAX(val, 0);
388
0
        val = L_MIN(val, 255);
389
0
        numaAddNumber(na, val);
390
0
    }
391
0
    for (i = maxval + 1; i < 256; i++)
392
0
        numaAddNumber(na, 255);
393
394
0
    return na;
395
0
}
396
397
398
/*-------------------------------------------------------------*
399
 *                      Contrast enhancement                   *
400
 *-------------------------------------------------------------*/
401
/*!
402
 * \brief   pixContrastTRC()
403
 *
404
 * \param[in]    pixd     [optional] null or equal to pixs
405
 * \param[in]    pixs     8 or 32 bpp; or 2, 4 or 8 bpp with colormap
406
 * \param[in]    factor   0.0 is no enhancement
407
 * \return  pixd always
408
 *
409
 * <pre>
410
 * Notes:
411
 *      (1) pixd must either be null or equal to pixs.
412
 *          For in-place operation, set pixd == pixs:
413
 *             pixContrastTRC(pixs, pixs, ...);
414
 *          To get a new image, set pixd == null:
415
 *             pixd = pixContrastTRC(NULL, pixs, ...);
416
 *      (2) If pixs is colormapped, the colormap is transformed,
417
 *          either in-place or in a copy of pixs.
418
 *      (3) Contrast is enhanced by mapping each color component
419
 *          using an atan function with maximum slope at 127.
420
 *          Pixels below 127 are lowered in intensity and pixels
421
 *          above 127 are increased.
422
 *      (4) The useful range for the contrast factor is scaled to
423
 *          be in (0.0 to 1.0), but larger values can also be used.
424
 *      (5) If factor == 0.0, no enhancement is performed; return a copy
425
 *          unless in-place, in which case this is a no-op.
426
 *      (6) For color images that are not colormapped, the mapping
427
 *          is applied to each component.
428
 * </pre>
429
 */
430
PIX *
431
pixContrastTRC(PIX       *pixd,
432
               PIX       *pixs,
433
               l_float32  factor)
434
0
{
435
0
l_int32   d;
436
0
NUMA     *nac;
437
0
PIXCMAP  *cmap;
438
439
0
    if (!pixs)
440
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
441
0
    if (pixd && (pixd != pixs))
442
0
        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
443
0
    if (factor < 0.0) {
444
0
        L_WARNING("factor must be >= 0.0; using 0.0\n", __func__);
445
0
        factor = 0.0;
446
0
    }
447
0
    if (factor == 0.0)
448
0
        return pixCopy(pixd, pixs);
449
450
0
    cmap = pixGetColormap(pixs);
451
0
    d = pixGetDepth(pixs);
452
0
    if (!cmap && d != 8 && d != 32)
453
0
        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, pixd);
454
455
0
    if (!pixd)  /* start with a copy if not in-place */
456
0
        pixd = pixCopy(NULL, pixs);
457
458
0
    if (cmap) {
459
0
        pixcmapContrastTRC(pixGetColormap(pixd), factor);
460
0
        return pixd;
461
0
    }
462
463
        /* pixd is 8 or 32 bpp */
464
0
    if ((nac = numaContrastTRC(factor)) == NULL)
465
0
        return (PIX *)ERROR_PTR("nac not made", __func__, pixd);
466
0
    pixTRCMap(pixd, NULL, nac);
467
0
    numaDestroy(&nac);
468
469
0
    return pixd;
470
0
}
471
472
473
/*!
474
 * \brief   pixContrastTRCMasked()
475
 *
476
 * \param[in]    pixd     [optional] null or equal to pixs
477
 * \param[in]    pixs     8 or 32 bpp; or 2, 4 or 8 bpp with colormap
478
 * \param[in]    pixm     [optional] null or 1 bpp
479
 * \param[in]    factor   0.0 is no enhancement
480
 * \return  pixd always
481
 *
482
 * <pre>
483
 * Notes:
484
 *      (1) Same as pixContrastTRC() except mapping is optionally over
485
 *          a subset of pixels described by pixm.
486
 *      (2) Masking does not work for colormapped images.
487
 *      (3) See pixContrastTRC() for details on how to use the parameters.
488
 * </pre>
489
 */
490
PIX *
491
pixContrastTRCMasked(PIX       *pixd,
492
                     PIX       *pixs,
493
                     PIX       *pixm,
494
                     l_float32  factor)
495
0
{
496
0
l_int32  d;
497
0
NUMA    *nac;
498
499
0
    if (!pixm)
500
0
        return pixContrastTRC(pixd, pixs, factor);
501
502
0
    if (!pixs)
503
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
504
0
    if (pixGetColormap(pixs))
505
0
        return (PIX *)ERROR_PTR("invalid: pixs has a colormap", __func__, pixd);
506
0
    if (pixd && (pixd != pixs))
507
0
        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
508
0
    d = pixGetDepth(pixs);
509
0
    if (d != 8 && d != 32)
510
0
        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, pixd);
511
512
0
    if (factor < 0.0) {
513
0
        L_WARNING("factor must be >= 0.0; using 0.0\n", __func__);
514
0
        factor = 0.0;
515
0
    }
516
0
    if (factor == 0.0)
517
0
        return pixCopy(pixd, pixs);
518
519
0
    if (!pixd)  /* start with a copy if not in-place */
520
0
        pixd = pixCopy(NULL, pixs);
521
522
0
    if ((nac = numaContrastTRC(factor)) == NULL)
523
0
        return (PIX *)ERROR_PTR("nac not made", __func__, pixd);
524
0
    pixTRCMap(pixd, pixm, nac);
525
0
    numaDestroy(&nac);
526
527
0
    return pixd;
528
0
}
529
530
531
/*!
532
 * \brief   numaContrastTRC()
533
 *
534
 * \param[in]    factor   generally between 0.0 [no enhancement]
535
 *                        and 1.0, but can be larger than 1.0
536
 * \return  na, or NULL on error
537
 *
538
 * <pre>
539
 * Notes:
540
 *      (1) The mapping is monotonic increasing, where 0 is mapped
541
 *          to 0 and 255 is mapped to 255.
542
 *      (2) As 'factor' is increased from 0.0 (where the mapping is linear),
543
 *          the map gets closer to its limit as a step function that
544
 *          jumps from 0 to 255 at the center (input value = 127).
545
 * </pre>
546
 */
547
NUMA *
548
numaContrastTRC(l_float32  factor)
549
0
{
550
0
l_int32    i, val;
551
0
l_float64  x, ymax, ymin, dely, scale;
552
0
NUMA      *na;
553
554
0
    if (factor < 0.0) {
555
0
        L_WARNING("factor must be >= 0.0; using 0.0; no enhancement\n",
556
0
                  __func__);
557
0
        factor = 0.0;
558
0
    }
559
0
    if (factor == 0.0)
560
0
        return numaMakeSequence(0, 1, 256);  /* linear map */
561
562
0
    scale = EnhanceScaleFactor;
563
0
    ymax = atan((l_float64)(1.0 * factor * scale));
564
0
    ymin = atan((l_float64)(-127. * factor * scale / 128.));
565
0
    dely = ymax - ymin;
566
0
    na = numaCreate(256);
567
0
    for (i = 0; i < 256; i++) {
568
0
        x = (l_float64)i;
569
0
        val = (l_int32)((255. / dely) *
570
0
             (-ymin + atan((l_float64)(factor * scale * (x - 127.) / 128.))) +
571
0
                 0.5);
572
0
        numaAddNumber(na, val);
573
0
    }
574
575
0
    return na;
576
0
}
577
578
579
/*-------------------------------------------------------------*
580
 *                     Histogram equalization                  *
581
 *-------------------------------------------------------------*/
582
/*!
583
 * \brief   pixEqualizeTRC()
584
 *
585
 * \param[in]    pixd     [optional] null or equal to pixs
586
 * \param[in]    pixs     8 bpp gray, 32 bpp rgb, or colormapped
587
 * \param[in]    fract    fraction of equalization movement of pixel values
588
 * \param[in]    factor   subsampling factor; integer >= 1
589
 * \return  pixd, or NULL on error
590
 *
591
 * <pre>
592
 * Notes:
593
 *      (1) pixd must either be null or equal to pixs.
594
 *          For in-place operation, set pixd == pixs:
595
 *             pixEqualizeTRC(pixs, pixs, ...);
596
 *          To get a new image, set pixd == null:
597
 *             pixd = pixEqualizeTRC(NULL, pixs, ...);
598
 *      (2) In histogram equalization, a tone reproduction curve
599
 *          mapping is used to make the number of pixels at each
600
 *          intensity equal.
601
 *      (3) If fract == 0.0, no equalization is performed; return a copy
602
 *          unless in-place, in which case this is a no-op.
603
 *          If fract == 1.0, equalization is complete.
604
 *      (4) Set the subsampling factor > 1 to reduce the amount of computation.
605
 *      (5) If pixs is colormapped, the colormap is removed and
606
 *          converted to rgb or grayscale.
607
 *      (6) If pixs has color, equalization is done in each channel
608
 *          separately.
609
 *      (7) Note that even if there is a colormap, we can get an
610
 *          in-place operation because the intermediate image pixt
611
 *          is copied back to pixs (which for in-place is the same
612
 *          as pixd).
613
 * </pre>
614
 */
615
PIX *
616
pixEqualizeTRC(PIX       *pixd,
617
               PIX       *pixs,
618
               l_float32  fract,
619
               l_int32    factor)
620
0
{
621
0
l_int32   d;
622
0
NUMA     *na;
623
0
PIX      *pixt, *pix8;
624
0
PIXCMAP  *cmap;
625
626
0
    if (!pixs)
627
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
628
0
    if (pixd && (pixd != pixs))
629
0
        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
630
0
    cmap = pixGetColormap(pixs);
631
0
    d = pixGetDepth(pixs);
632
0
    if (d != 8 && d != 32 && !cmap)
633
0
        return (PIX *)ERROR_PTR("pixs not 8/32 bpp or cmapped", __func__, NULL);
634
0
    if (fract < 0.0 || fract > 1.0)
635
0
        return (PIX *)ERROR_PTR("fract not in [0.0 ... 1.0]", __func__, NULL);
636
0
    if (factor < 1)
637
0
        return (PIX *)ERROR_PTR("sampling factor < 1", __func__, NULL);
638
639
0
    if (fract == 0.0)
640
0
        return pixCopy(pixd, pixs);
641
642
        /* If there is a colormap, remove it. */
643
0
    if (cmap)
644
0
        pixt = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
645
0
    else
646
0
        pixt = pixClone(pixs);
647
648
        /* Make a copy if necessary */
649
0
    pixd = pixCopy(pixd, pixt);
650
0
    pixDestroy(&pixt);
651
652
0
    d = pixGetDepth(pixd);
653
0
    if (d == 8) {
654
0
        na = numaEqualizeTRC(pixd, fract, factor);
655
0
        pixTRCMap(pixd, NULL, na);
656
0
        numaDestroy(&na);
657
0
    } else {  /* 32 bpp */
658
0
        pix8 = pixGetRGBComponent(pixd, COLOR_RED);
659
0
        na = numaEqualizeTRC(pix8, fract, factor);
660
0
        pixTRCMap(pix8, NULL, na);
661
0
        pixSetRGBComponent(pixd, pix8, COLOR_RED);
662
0
        numaDestroy(&na);
663
0
        pixDestroy(&pix8);
664
0
        pix8 = pixGetRGBComponent(pixd, COLOR_GREEN);
665
0
        na = numaEqualizeTRC(pix8, fract, factor);
666
0
        pixTRCMap(pix8, NULL, na);
667
0
        pixSetRGBComponent(pixd, pix8, COLOR_GREEN);
668
0
        numaDestroy(&na);
669
0
        pixDestroy(&pix8);
670
0
        pix8 = pixGetRGBComponent(pixd, COLOR_BLUE);
671
0
        na = numaEqualizeTRC(pix8, fract, factor);
672
0
        pixTRCMap(pix8, NULL, na);
673
0
        pixSetRGBComponent(pixd, pix8, COLOR_BLUE);
674
0
        numaDestroy(&na);
675
0
        pixDestroy(&pix8);
676
0
    }
677
678
0
    return pixd;
679
0
}
680
681
682
/*!
683
 * \brief   numaEqualizeTRC()
684
 *
685
 * \param[in]    pix     8 bpp, no colormap
686
 * \param[in]    fract   fraction of equalization movement of pixel values
687
 * \param[in]    factor  subsampling factor; integer >= 1
688
 * \return  nad, or NULL on error
689
 *
690
 * <pre>
691
 * Notes:
692
 *      (1) If fract == 0.0, no equalization will be performed.
693
 *          If fract == 1.0, equalization is complete.
694
 *      (2) Set the subsampling factor > 1 to reduce the amount of computation.
695
 *      (3) The map is returned as a numa with 256 values, specifying
696
 *          the equalized value (array value) for every input value
697
 *          (the array index).
698
 * </pre>
699
 */
700
NUMA *
701
numaEqualizeTRC(PIX       *pix,
702
                l_float32  fract,
703
                l_int32    factor)
704
0
{
705
0
l_int32    iin, iout, itarg;
706
0
l_float32  val, sum;
707
0
NUMA      *nah, *nasum, *nad;
708
709
0
    if (!pix)
710
0
        return (NUMA *)ERROR_PTR("pix not defined", __func__, NULL);
711
0
    if (pixGetDepth(pix) != 8)
712
0
        return (NUMA *)ERROR_PTR("pix not 8 bpp", __func__, NULL);
713
0
    if (fract < 0.0 || fract > 1.0)
714
0
        return (NUMA *)ERROR_PTR("fract not in [0.0 ... 1.0]", __func__, NULL);
715
0
    if (factor < 1)
716
0
        return (NUMA *)ERROR_PTR("sampling factor < 1", __func__, NULL);
717
718
0
    if (fract == 0.0)
719
0
        L_WARNING("fract = 0.0; no equalization requested\n", __func__);
720
721
0
    if ((nah = pixGetGrayHistogram(pix, factor)) == NULL)
722
0
        return (NUMA *)ERROR_PTR("histogram not made", __func__, NULL);
723
0
    numaGetSum(nah, &sum);
724
0
    nasum = numaGetPartialSums(nah);
725
726
0
    nad = numaCreate(256);
727
0
    for (iin = 0; iin < 256; iin++) {
728
0
        numaGetFValue(nasum, iin, &val);
729
0
        itarg = (l_int32)(255. * val / sum + 0.5);
730
0
        iout = iin + (l_int32)(fract * (itarg - iin));
731
0
        iout = L_MIN(iout, 255);  /* to be safe */
732
0
        numaAddNumber(nad, iout);
733
0
    }
734
735
0
    numaDestroy(&nah);
736
0
    numaDestroy(&nasum);
737
0
    return nad;
738
0
}
739
740
741
/*-------------------------------------------------------------*
742
 *                       Generic TRC mapping                   *
743
 *-------------------------------------------------------------*/
744
/*!
745
 * \brief   pixTRCMap()
746
 *
747
 * \param[in]    pixs    8 grayscale or 32 bpp rgb; not colormapped
748
 * \param[in]    pixm    [optional] 1 bpp mask
749
 * \param[in]    na      mapping array
750
 * \return  0 if OK, 1 on error
751
 *
752
 * <pre>
753
 * Notes:
754
 *      (1) This operation is in-place on pixs.
755
 *      (2) For 32 bpp, this applies the same map to each of the r,g,b
756
 *          components.
757
 *      (3) The mapping array is of size 256, and it maps the input
758
 *          index into values in the range [0, 255].
759
 *      (4) If defined, the optional 1 bpp mask pixm has its origin
760
 *          aligned with pixs, and the map function is applied only
761
 *          to pixels in pixs under the fg of pixm.
762
 *      (5) For 32 bpp, this does not save the alpha channel.
763
 * </pre>
764
 */
765
l_int32
766
pixTRCMap(PIX   *pixs,
767
          PIX   *pixm,
768
          NUMA  *na)
769
0
{
770
0
l_int32    w, h, d, wm, hm, wpl, wplm, i, j, sval8, dval8;
771
0
l_uint32   sval32, dval32;
772
0
l_uint32  *data, *datam, *line, *linem, *tab;
773
774
0
    if (!pixs)
775
0
        return ERROR_INT("pixs not defined", __func__, 1);
776
0
    if (pixGetColormap(pixs))
777
0
        return ERROR_INT("pixs is colormapped", __func__, 1);
778
0
    if (!na)
779
0
        return ERROR_INT("na not defined", __func__, 1);
780
0
    if (numaGetCount(na) != 256)
781
0
        return ERROR_INT("na not of size 256", __func__, 1);
782
0
    pixGetDimensions(pixs, &w, &h, &d);
783
0
    if (d != 8 && d != 32)
784
0
        return ERROR_INT("pixs not 8 or 32 bpp", __func__, 1);
785
0
    if (pixm) {
786
0
        if (pixGetDepth(pixm) != 1)
787
0
            return ERROR_INT("pixm not 1 bpp", __func__, 1);
788
0
    }
789
790
0
    tab = (l_uint32 *)numaGetIArray(na);  /* get the array for efficiency */
791
0
    wpl = pixGetWpl(pixs);
792
0
    data = pixGetData(pixs);
793
0
    if (!pixm) {
794
0
        if (d == 8) {
795
0
            for (i = 0; i < h; i++) {
796
0
                line = data + i * wpl;
797
0
                for (j = 0; j < w; j++) {
798
0
                    sval8 = GET_DATA_BYTE(line, j);
799
0
                    dval8 = tab[sval8];
800
0
                    SET_DATA_BYTE(line, j, dval8);
801
0
                }
802
0
            }
803
0
        } else {  /* d == 32 */
804
0
            for (i = 0; i < h; i++) {
805
0
                line = data + i * wpl;
806
0
                for (j = 0; j < w; j++) {
807
0
                    sval32 = *(line + j);
808
0
                    dval32 =
809
0
                        tab[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
810
0
                        tab[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
811
0
                        tab[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
812
0
                    *(line + j) = dval32;
813
0
                }
814
0
            }
815
0
        }
816
0
    } else {
817
0
        datam = pixGetData(pixm);
818
0
        wplm = pixGetWpl(pixm);
819
0
        pixGetDimensions(pixm, &wm, &hm, NULL);
820
0
        if (d == 8) {
821
0
            for (i = 0; i < h; i++) {
822
0
                if (i >= hm)
823
0
                    break;
824
0
                line = data + i * wpl;
825
0
                linem = datam + i * wplm;
826
0
                for (j = 0; j < w; j++) {
827
0
                    if (j >= wm)
828
0
                        break;
829
0
                    if (GET_DATA_BIT(linem, j) == 0)
830
0
                        continue;
831
0
                    sval8 = GET_DATA_BYTE(line, j);
832
0
                    dval8 = tab[sval8];
833
0
                    SET_DATA_BYTE(line, j, dval8);
834
0
                }
835
0
            }
836
0
        } else {  /* d == 32 */
837
0
            for (i = 0; i < h; i++) {
838
0
                if (i >= hm)
839
0
                    break;
840
0
                line = data + i * wpl;
841
0
                linem = datam + i * wplm;
842
0
                for (j = 0; j < w; j++) {
843
0
                    if (j >= wm)
844
0
                        break;
845
0
                    if (GET_DATA_BIT(linem, j) == 0)
846
0
                        continue;
847
0
                    sval32 = *(line + j);
848
0
                    dval32 =
849
0
                        tab[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
850
0
                        tab[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
851
0
                        tab[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
852
0
                    *(line + j) = dval32;
853
0
                }
854
0
            }
855
0
        }
856
0
    }
857
858
0
    LEPT_FREE(tab);
859
0
    return 0;
860
0
}
861
862
863
/*!
864
 * \brief   pixTRCMapGeneral()
865
 *
866
 * \param[in]    pixs             32 bpp rgb; not colormapped
867
 * \param[in]    pixm             [optional] 1 bpp mask
868
 * \param[in]    nar, nag, nab    mapping arrays
869
 * \return  0 if OK, 1 on error
870
 *
871
 * <pre>
872
 * Notes:
873
 *      (1) This operation is in-place on %pixs.
874
 *      (2) Each of the r,g,b mapping arrays is of size 256. They map the
875
 *          input value for that color component into values in the
876
 *          range [0, 255].
877
 *      (3) In the special case where the r, g and b mapping arrays are
878
 *          all the same, call pixTRCMap() instead.
879
 *      (4) If defined, the optional 1 bpp mask %pixm has its origin
880
 *          aligned with %pixs, and the map function is applied only
881
 *          to pixels in %pixs under the fg of pixm.
882
 *      (5) The alpha channel is not saved.
883
 * </pre>
884
 */
885
l_int32
886
pixTRCMapGeneral(PIX   *pixs,
887
                 PIX   *pixm,
888
                 NUMA  *nar,
889
                 NUMA  *nag,
890
                 NUMA  *nab)
891
0
{
892
0
l_int32    w, h, wm, hm, wpl, wplm, i, j;
893
0
l_uint32   sval32, dval32;
894
0
l_uint32  *data, *datam, *line, *linem, *tabr, *tabg, *tabb;
895
896
0
    if (!pixs || pixGetDepth(pixs) != 32)
897
0
        return ERROR_INT("pixs not defined or not 32 bpp", __func__, 1);
898
0
    if (pixm && pixGetDepth(pixm) != 1)
899
0
        return ERROR_INT("pixm defined and not 1 bpp", __func__, 1);
900
0
    if (!nar || !nag || !nab)
901
0
        return ERROR_INT("na{r,g,b} not all defined", __func__, 1);
902
0
    if (numaGetCount(nar) != 256 || numaGetCount(nag) != 256 ||
903
0
        numaGetCount(nab) != 256)
904
0
        return ERROR_INT("na{r,g,b} not all of size 256", __func__, 1);
905
906
        /* Get the arrays for efficiency */
907
0
    tabr = (l_uint32 *)numaGetIArray(nar);
908
0
    tabg = (l_uint32 *)numaGetIArray(nag);
909
0
    tabb = (l_uint32 *)numaGetIArray(nab);
910
0
    pixGetDimensions(pixs, &w, &h, NULL);
911
0
    wpl = pixGetWpl(pixs);
912
0
    data = pixGetData(pixs);
913
0
    if (!pixm) {
914
0
        for (i = 0; i < h; i++) {
915
0
            line = data + i * wpl;
916
0
            for (j = 0; j < w; j++) {
917
0
                sval32 = *(line + j);
918
0
                dval32 =
919
0
                    tabr[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
920
0
                    tabg[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
921
0
                    tabb[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
922
0
                *(line + j) = dval32;
923
0
            }
924
0
        }
925
0
    } else {
926
0
        datam = pixGetData(pixm);
927
0
        wplm = pixGetWpl(pixm);
928
0
        pixGetDimensions(pixm, &wm, &hm, NULL);
929
0
        for (i = 0; i < h; i++) {
930
0
            if (i >= hm)
931
0
                break;
932
0
            line = data + i * wpl;
933
0
            linem = datam + i * wplm;
934
0
            for (j = 0; j < w; j++) {
935
0
                if (j >= wm)
936
0
                    break;
937
0
                if (GET_DATA_BIT(linem, j) == 0)
938
0
                    continue;
939
0
                sval32 = *(line + j);
940
0
                dval32 =
941
0
                    tabr[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
942
0
                    tabg[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
943
0
                    tabb[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
944
0
                *(line + j) = dval32;
945
0
            }
946
0
        }
947
0
    }
948
949
0
    LEPT_FREE(tabr);
950
0
    LEPT_FREE(tabg);
951
0
    LEPT_FREE(tabb);
952
0
    return 0;
953
0
}
954
955
956
957
/*-----------------------------------------------------------------------*
958
 *                             Unsharp masking                           *
959
 *-----------------------------------------------------------------------*/
960
/*!
961
 * \brief   pixUnsharpMasking()
962
 *
963
 * \param[in]    pixs       all depths except 1 bpp; with or without colormaps
964
 * \param[in]    halfwidth  "half-width" of smoothing filter
965
 * \param[in]    fract      fraction of edge added back into image
966
 * \return  pixd, or NULL on error
967
 *
968
 * <pre>
969
 * Notes:
970
 *      (1) We use symmetric smoothing filters of odd dimension,
971
 *          typically use sizes of 3, 5, 7, etc.  The %halfwidth parameter
972
 *          for these is (size - 1)/2; i.e., 1, 2, 3, etc.
973
 *      (2) The fract parameter is typically taken in the
974
 *          range:  0.2 < fract < 0.7
975
 *      (3) Returns a clone if no sharpening is requested.
976
 * </pre>
977
 */
978
PIX *
979
pixUnsharpMasking(PIX       *pixs,
980
                  l_int32    halfwidth,
981
                  l_float32  fract)
982
0
{
983
0
l_int32  d;
984
0
PIX     *pix1, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;
985
986
0
    if (!pixs || (pixGetDepth(pixs) == 1))
987
0
        return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", __func__, NULL);
988
0
    if (fract <= 0.0 || halfwidth <= 0) {
989
0
        L_WARNING("no sharpening requested; clone returned\n", __func__);
990
0
        return pixClone(pixs);
991
0
    }
992
993
0
    if (halfwidth == 1 || halfwidth == 2)
994
0
        return pixUnsharpMaskingFast(pixs, halfwidth, fract, L_BOTH_DIRECTIONS);
995
996
        /* Remove colormap; clone if possible; result is either 8 or 32 bpp */
997
0
    if ((pix1 = pixConvertTo8Or32(pixs, L_CLONE, 0)) == NULL)
998
0
        return (PIX *)ERROR_PTR("pix1 not made", __func__, NULL);
999
1000
        /* Sharpen */
1001
0
    d = pixGetDepth(pix1);
1002
0
    if (d == 8) {
1003
0
        pixd = pixUnsharpMaskingGray(pix1, halfwidth, fract);
1004
0
    } else {  /* d == 32 */
1005
0
        pixr = pixGetRGBComponent(pix1, COLOR_RED);
1006
0
        pixrs = pixUnsharpMaskingGray(pixr, halfwidth, fract);
1007
0
        pixDestroy(&pixr);
1008
0
        pixg = pixGetRGBComponent(pix1, COLOR_GREEN);
1009
0
        pixgs = pixUnsharpMaskingGray(pixg, halfwidth, fract);
1010
0
        pixDestroy(&pixg);
1011
0
        pixb = pixGetRGBComponent(pix1, COLOR_BLUE);
1012
0
        pixbs = pixUnsharpMaskingGray(pixb, halfwidth, fract);
1013
0
        pixDestroy(&pixb);
1014
0
        pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
1015
0
        pixDestroy(&pixrs);
1016
0
        pixDestroy(&pixgs);
1017
0
        pixDestroy(&pixbs);
1018
0
        if (pixGetSpp(pixs) == 4)
1019
0
            pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
1020
0
    }
1021
1022
0
    pixDestroy(&pix1);
1023
0
    return pixd;
1024
0
}
1025
1026
1027
/*!
1028
 * \brief   pixUnsharpMaskingGray()
1029
 *
1030
 * \param[in]    pixs       8 bpp; no colormap
1031
 * \param[in]    halfwidth  "half-width" of smoothing filter
1032
 * \param[in]    fract      fraction of edge added back into image
1033
 * \return  pixd, or NULL on error
1034
 *
1035
 * <pre>
1036
 * Notes:
1037
 *      (1) We use symmetric smoothing filters of odd dimension,
1038
 *          typically use sizes of 3, 5, 7, etc.  The %halfwidth parameter
1039
 *          for these is (size - 1)/2; i.e., 1, 2, 3, etc.
1040
 *      (2) The fract parameter is typically taken in the range:
1041
 *          0.2 < fract < 0.7
1042
 *      (3) Returns a clone if no sharpening is requested.
1043
 * </pre>
1044
 */
1045
PIX *
1046
pixUnsharpMaskingGray(PIX       *pixs,
1047
                      l_int32    halfwidth,
1048
                      l_float32  fract)
1049
0
{
1050
0
l_int32  w, h, d;
1051
0
PIX     *pixc, *pixd;
1052
0
PIXACC  *pixacc;
1053
1054
0
    if (!pixs)
1055
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1056
0
    pixGetDimensions(pixs, &w, &h, &d);
1057
0
    if (d != 8 || pixGetColormap(pixs) != NULL)
1058
0
        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", __func__, NULL);
1059
0
    if (fract <= 0.0 || halfwidth <= 0) {
1060
0
        L_WARNING("no sharpening requested; clone returned\n", __func__);
1061
0
        return pixClone(pixs);
1062
0
    }
1063
0
    if (halfwidth == 1 || halfwidth == 2)
1064
0
        return pixUnsharpMaskingGrayFast(pixs, halfwidth, fract,
1065
0
                                         L_BOTH_DIRECTIONS);
1066
1067
0
    if ((pixc = pixBlockconvGray(pixs, NULL, halfwidth, halfwidth)) == NULL)
1068
0
        return (PIX *)ERROR_PTR("pixc not made", __func__, NULL);
1069
1070
        /* Steps:
1071
         *    (1) edge image is pixs - pixc  (this is highpass part)
1072
         *    (2) multiply edge image by fract
1073
         *    (3) add fraction of edge to pixs
1074
         *
1075
         * To show how this is done with both interfaces to arithmetic
1076
         * on integer Pix, here is the implementation in the lower-level
1077
         * function calls:
1078
         *    pixt = pixInitAccumulate(w, h, 0x10000000)) == NULL)
1079
         *    pixAccumulate(pixt, pixs, L_ARITH_ADD);
1080
         *    pixAccumulate(pixt, pixc, L_ARITH_SUBTRACT);
1081
         *    pixMultConstAccumulate(pixt, fract, 0x10000000);
1082
         *    pixAccumulate(pixt, pixs, L_ARITH_ADD);
1083
         *    pixd = pixFinalAccumulate(pixt, 0x10000000, 8)) == NULL)
1084
         *    pixDestroy(&pixt);
1085
         *
1086
         * The code below does the same thing using the Pixacc accumulator,
1087
         * hiding the details of the offset that is needed for subtraction.
1088
         */
1089
0
    pixacc = pixaccCreate(w, h, 1);
1090
0
    pixaccAdd(pixacc, pixs);
1091
0
    pixaccSubtract(pixacc, pixc);
1092
0
    pixaccMultConst(pixacc, fract);
1093
0
    pixaccAdd(pixacc, pixs);
1094
0
    pixd = pixaccFinal(pixacc, 8);
1095
0
    pixaccDestroy(&pixacc);
1096
1097
0
    pixDestroy(&pixc);
1098
0
    return pixd;
1099
0
}
1100
1101
1102
/*!
1103
 * \brief   pixUnsharpMaskingFast()
1104
 *
1105
 * \param[in]    pixs       all depths except 1 bpp; with or without colormaps
1106
 * \param[in]    halfwidth  "half-width" of smoothing filter; 1 and 2 only
1107
 * \param[in]    fract      fraction of high frequency added to image
1108
 * \param[in]    direction  L_HORIZ, L_VERT, L_BOTH_DIRECTIONS
1109
 * \return  pixd, or NULL on error
1110
 *
1111
 * <pre>
1112
 * Notes:
1113
 *      (1) The fast version uses separable 1-D filters directly on
1114
 *          the input image.  The halfwidth is either 1 (full width = 3)
1115
 *          or 2 (full width = 5).
1116
 *      (2) The fract parameter is typically taken in the
1117
 *            range:  0.2 < fract < 0.7
1118
 *      (3) To skip horizontal sharpening, use %fracth = 0.0; ditto for %fractv
1119
 *      (4) For one dimensional filtering (as an example):
1120
 *          For %halfwidth = 1, the low-pass filter is
1121
 *              L:    1/3    1/3   1/3
1122
 *          and the high-pass filter is
1123
 *              H = I - L:   -1/3   2/3   -1/3
1124
 *          For %halfwidth = 2, the low-pass filter is
1125
 *              L:    1/5    1/5   1/5    1/5    1/5
1126
 *          and the high-pass filter is
1127
 *              H = I - L:   -1/5  -1/5   4/5  -1/5   -1/5
1128
 *          The new sharpened pixel value is found by adding some fraction
1129
 *          of the high-pass filter value (which sums to 0) to the
1130
 *          initial pixel value:
1131
 *              N = I + fract * H
1132
 *      (5) For 2D, the sharpening filter is not separable, because the
1133
 *          vertical filter depends on the horizontal location relative
1134
 *          to the filter origin, and v.v.   So we either do the full
1135
 *          2D filter (for %halfwidth == 1) or do the low-pass
1136
 *          convolution separably and then compose with the original pix.
1137
 *      (6) Returns a clone if no sharpening is requested.
1138
 * </pre>
1139
 */
1140
PIX *
1141
pixUnsharpMaskingFast(PIX       *pixs,
1142
                      l_int32    halfwidth,
1143
                      l_float32  fract,
1144
                      l_int32    direction)
1145
0
{
1146
0
l_int32  d;
1147
0
PIX     *pixt, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;
1148
1149
0
    if (!pixs || (pixGetDepth(pixs) == 1))
1150
0
        return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", __func__, NULL);
1151
0
    if (fract <= 0.0 || halfwidth <= 0) {
1152
0
        L_WARNING("no sharpening requested; clone returned\n", __func__);
1153
0
        return pixClone(pixs);
1154
0
    }
1155
0
    if (halfwidth != 1 && halfwidth != 2)
1156
0
        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", __func__, NULL);
1157
0
    if (direction != L_HORIZ && direction != L_VERT &&
1158
0
        direction != L_BOTH_DIRECTIONS)
1159
0
        return (PIX *)ERROR_PTR("invalid direction", __func__, NULL);
1160
1161
        /* Remove colormap; clone if possible; result is either 8 or 32 bpp */
1162
0
    if ((pixt = pixConvertTo8Or32(pixs, L_CLONE, 0)) == NULL)
1163
0
        return (PIX *)ERROR_PTR("pixt not made", __func__, NULL);
1164
1165
        /* Sharpen */
1166
0
    d = pixGetDepth(pixt);
1167
0
    if (d == 8) {
1168
0
        pixd = pixUnsharpMaskingGrayFast(pixt, halfwidth, fract, direction);
1169
0
    } else {  /* d == 32 */
1170
0
        pixr = pixGetRGBComponent(pixs, COLOR_RED);
1171
0
        pixrs = pixUnsharpMaskingGrayFast(pixr, halfwidth, fract, direction);
1172
0
        pixDestroy(&pixr);
1173
0
        pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
1174
0
        pixgs = pixUnsharpMaskingGrayFast(pixg, halfwidth, fract, direction);
1175
0
        pixDestroy(&pixg);
1176
0
        pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
1177
0
        pixbs = pixUnsharpMaskingGrayFast(pixb, halfwidth, fract, direction);
1178
0
        pixDestroy(&pixb);
1179
0
        pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
1180
0
        if (pixGetSpp(pixs) == 4)
1181
0
            pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
1182
0
        pixDestroy(&pixrs);
1183
0
        pixDestroy(&pixgs);
1184
0
        pixDestroy(&pixbs);
1185
0
    }
1186
1187
0
    pixDestroy(&pixt);
1188
0
    return pixd;
1189
0
}
1190
1191
1192
1193
/*!
1194
 * \brief   pixUnsharpMaskingGrayFast()
1195
 *
1196
 * \param[in]    pixs       8 bpp; no colormap
1197
 * \param[in]    halfwidth  "half-width" of smoothing filter: 1 or 2
1198
 * \param[in]    fract      fraction of high frequency added to image
1199
 * \param[in]    direction  L_HORIZ, L_VERT, L_BOTH_DIRECTIONS
1200
 * \return  pixd, or NULL on error
1201
 *
1202
 * <pre>
1203
 * Notes:
1204
 *      (1) For usage and explanation of the algorithm, see notes
1205
 *          in pixUnsharpMaskingFast().
1206
 *      (2) Returns a clone if no sharpening is requested.
1207
 * </pre>
1208
 */
1209
PIX *
1210
pixUnsharpMaskingGrayFast(PIX       *pixs,
1211
                          l_int32    halfwidth,
1212
                          l_float32  fract,
1213
                          l_int32    direction)
1214
0
{
1215
0
PIX  *pixd;
1216
1217
0
    if (!pixs)
1218
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1219
0
    if (pixGetDepth(pixs) != 8 || pixGetColormap(pixs) != NULL)
1220
0
        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", __func__, NULL);
1221
0
    if (fract <= 0.0 || halfwidth <= 0) {
1222
0
        L_WARNING("no sharpening requested; clone returned\n", __func__);
1223
0
        return pixClone(pixs);
1224
0
    }
1225
0
    if (halfwidth != 1 && halfwidth != 2)
1226
0
        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", __func__, NULL);
1227
0
    if (direction != L_HORIZ && direction != L_VERT &&
1228
0
        direction != L_BOTH_DIRECTIONS)
1229
0
        return (PIX *)ERROR_PTR("invalid direction", __func__, NULL);
1230
1231
0
    if (direction != L_BOTH_DIRECTIONS)
1232
0
        pixd = pixUnsharpMaskingGray1D(pixs, halfwidth, fract, direction);
1233
0
    else  /* 2D sharpening */
1234
0
        pixd = pixUnsharpMaskingGray2D(pixs, halfwidth, fract);
1235
1236
0
    return pixd;
1237
0
}
1238
1239
1240
/*!
1241
 * \brief   pixUnsharpMaskingGray1D()
1242
 *
1243
 * \param[in]    pixs        8 bpp; no colormap
1244
 * \param[in]    halfwidth   "half-width" of smoothing filter: 1 or 2
1245
 * \param[in]    fract       fraction of high frequency added to image
1246
 * \param[in]    direction   filtering direction; use L_HORIZ or L_VERT
1247
 * \return  pixd, or NULL on error
1248
 *
1249
 * <pre>
1250
 * Notes:
1251
 *      (1) For usage and explanation of the algorithm, see notes
1252
 *          in pixUnsharpMaskingFast().
1253
 *      (2) Returns a clone if no sharpening is requested.
1254
 * </pre>
1255
 */
1256
PIX *
1257
pixUnsharpMaskingGray1D(PIX       *pixs,
1258
                        l_int32    halfwidth,
1259
                        l_float32  fract,
1260
                        l_int32    direction)
1261
0
{
1262
0
l_int32    w, h, d, wpls, wpld, i, j, ival;
1263
0
l_uint32  *datas, *datad;
1264
0
l_uint32  *lines, *lines0, *lines1, *lines2, *lines3, *lines4, *lined;
1265
0
l_float32  val, a[5];
1266
0
PIX       *pixd;
1267
1268
0
    if (!pixs)
1269
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1270
0
    pixGetDimensions(pixs, &w, &h, &d);
1271
0
    if (d != 8 || pixGetColormap(pixs) != NULL)
1272
0
        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", __func__, NULL);
1273
0
    if (fract <= 0.0 || halfwidth <= 0) {
1274
0
        L_WARNING("no sharpening requested; clone returned\n", __func__);
1275
0
        return pixClone(pixs);
1276
0
    }
1277
0
    if (halfwidth != 1 && halfwidth != 2)
1278
0
        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", __func__, NULL);
1279
1280
        /* Initialize pixd with pixels from pixs that will not be
1281
         * set when computing the sharpened values. */
1282
0
    pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
1283
0
                         halfwidth, halfwidth);
1284
0
    datas = pixGetData(pixs);
1285
0
    datad = pixGetData(pixd);
1286
0
    wpls = pixGetWpl(pixs);
1287
0
    wpld = pixGetWpl(pixd);
1288
1289
0
    if (halfwidth == 1) {
1290
0
        a[0] = -fract / 3.0;
1291
0
        a[1] = 1.0 + fract * 2.0 / 3.0;
1292
0
        a[2] = a[0];
1293
0
    } else {  /* halfwidth == 2 */
1294
0
        a[0] = -fract / 5.0;
1295
0
        a[1] = a[0];
1296
0
        a[2] = 1.0 + fract * 4.0 / 5.0;
1297
0
        a[3] = a[0];
1298
0
        a[4] = a[0];
1299
0
    }
1300
1301
0
    if (direction == L_HORIZ) {
1302
0
        for (i = 0; i < h; i++) {
1303
0
            lines = datas + i * wpls;
1304
0
            lined = datad + i * wpld;
1305
0
            if (halfwidth == 1) {
1306
0
                for (j = 1; j < w - 1; j++) {
1307
0
                    val = a[0] * GET_DATA_BYTE(lines, j - 1) +
1308
0
                          a[1] * GET_DATA_BYTE(lines, j) +
1309
0
                          a[2] * GET_DATA_BYTE(lines, j + 1);
1310
0
                    ival = (l_int32)val;
1311
0
                    ival = L_MAX(0, ival);
1312
0
                    ival = L_MIN(255, ival);
1313
0
                    SET_DATA_BYTE(lined, j, ival);
1314
0
                }
1315
0
            } else {  /* halfwidth == 2 */
1316
0
                for (j = 2; j < w - 2; j++) {
1317
0
                    val = a[0] * GET_DATA_BYTE(lines, j - 2) +
1318
0
                          a[1] * GET_DATA_BYTE(lines, j - 1) +
1319
0
                          a[2] * GET_DATA_BYTE(lines, j) +
1320
0
                          a[3] * GET_DATA_BYTE(lines, j + 1) +
1321
0
                          a[4] * GET_DATA_BYTE(lines, j + 2);
1322
0
                    ival = (l_int32)val;
1323
0
                    ival = L_MAX(0, ival);
1324
0
                    ival = L_MIN(255, ival);
1325
0
                    SET_DATA_BYTE(lined, j, ival);
1326
0
                }
1327
0
            }
1328
0
        }
1329
0
    } else {  /* direction == L_VERT */
1330
0
        if (halfwidth == 1) {
1331
0
            for (i = 1; i < h - 1; i++) {
1332
0
                lines0 = datas + (i - 1) * wpls;
1333
0
                lines1 = datas + i * wpls;
1334
0
                lines2 = datas + (i + 1) * wpls;
1335
0
                lined = datad + i * wpld;
1336
0
                for (j = 0; j < w; j++) {
1337
0
                    val = a[0] * GET_DATA_BYTE(lines0, j) +
1338
0
                          a[1] * GET_DATA_BYTE(lines1, j) +
1339
0
                          a[2] * GET_DATA_BYTE(lines2, j);
1340
0
                    ival = (l_int32)val;
1341
0
                    ival = L_MAX(0, ival);
1342
0
                    ival = L_MIN(255, ival);
1343
0
                    SET_DATA_BYTE(lined, j, ival);
1344
0
                }
1345
0
            }
1346
0
        } else {  /* halfwidth == 2 */
1347
0
            for (i = 2; i < h - 2; i++) {
1348
0
                lines0 = datas + (i - 2) * wpls;
1349
0
                lines1 = datas + (i - 1) * wpls;
1350
0
                lines2 = datas + i * wpls;
1351
0
                lines3 = datas + (i + 1) * wpls;
1352
0
                lines4 = datas + (i + 2) * wpls;
1353
0
                lined = datad + i * wpld;
1354
0
                for (j = 0; j < w; j++) {
1355
0
                    val = a[0] * GET_DATA_BYTE(lines0, j) +
1356
0
                          a[1] * GET_DATA_BYTE(lines1, j) +
1357
0
                          a[2] * GET_DATA_BYTE(lines2, j) +
1358
0
                          a[3] * GET_DATA_BYTE(lines3, j) +
1359
0
                          a[4] * GET_DATA_BYTE(lines4, j);
1360
0
                    ival = (l_int32)val;
1361
0
                    ival = L_MAX(0, ival);
1362
0
                    ival = L_MIN(255, ival);
1363
0
                    SET_DATA_BYTE(lined, j, ival);
1364
0
                }
1365
0
            }
1366
0
        }
1367
0
    }
1368
1369
0
    return pixd;
1370
0
}
1371
1372
1373
/*!
1374
 * \brief   pixUnsharpMaskingGray2D()
1375
 *
1376
 * \param[in]    pixs       8 bpp; no colormap
1377
 * \param[in]    halfwidth  "half-width" of smoothing filter: 1 or 2
1378
 * \param[in]    fract      fraction of high frequency added to image
1379
 * \return  pixd, or NULL on error
1380
 *
1381
 * <pre>
1382
 * Notes:
1383
 *      (1) This is for %halfwidth == 1, 2.
1384
 *      (2) The lowpass filter is implemented separably.
1385
 *      (3) Returns a clone if no sharpening is requested.
1386
 * </pre>
1387
 */
1388
PIX *
1389
pixUnsharpMaskingGray2D(PIX       *pixs,
1390
                        l_int32    halfwidth,
1391
                        l_float32  fract)
1392
0
{
1393
0
l_int32     w, h, d, wpls, wpld, wplf, i, j, ival, sval;
1394
0
l_uint32   *datas, *datad, *lines, *lined;
1395
0
l_float32   val, norm;
1396
0
l_float32  *dataf, *linef, *linef0, *linef1, *linef2, *linef3, *linef4;
1397
0
PIX        *pixd;
1398
0
FPIX       *fpix;
1399
1400
0
    if (!pixs)
1401
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1402
0
    pixGetDimensions(pixs, &w, &h, &d);
1403
0
    if (d != 8 || pixGetColormap(pixs) != NULL)
1404
0
        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", __func__, NULL);
1405
0
    if (fract <= 0.0 || halfwidth <= 0) {
1406
0
        L_WARNING("no sharpening requested; clone returned\n", __func__);
1407
0
        return pixClone(pixs);
1408
0
    }
1409
0
    if (halfwidth != 1 && halfwidth != 2)
1410
0
        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", __func__, NULL);
1411
1412
0
    if ((pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
1413
0
                              halfwidth, halfwidth)) == NULL)
1414
0
        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
1415
0
    datad = pixGetData(pixd);
1416
0
    wpld = pixGetWpl(pixd);
1417
0
    datas = pixGetData(pixs);
1418
0
    wpls = pixGetWpl(pixs);
1419
1420
        /* Do the low pass separably.  Store the result of horizontal
1421
         * smoothing in an intermediate fpix.  */
1422
0
    if ((fpix = fpixCreate(w, h)) == NULL) {
1423
0
        pixDestroy(&pixd);
1424
0
        return (PIX *)ERROR_PTR("fpix not made", __func__, NULL);
1425
0
    }
1426
0
    dataf = fpixGetData(fpix);
1427
0
    wplf = fpixGetWpl(fpix);
1428
0
    if (halfwidth == 1) {
1429
0
        for (i = 0; i < h; i++) {
1430
0
            lines = datas + i * wpls;
1431
0
            linef = dataf + i * wplf;
1432
0
            for (j = 1; j < w - 1; j++) {
1433
0
                val = GET_DATA_BYTE(lines, j - 1) +
1434
0
                      GET_DATA_BYTE(lines, j) +
1435
0
                      GET_DATA_BYTE(lines, j + 1);
1436
0
                linef[j] = val;
1437
0
            }
1438
0
        }
1439
0
    } else {
1440
0
        for (i = 0; i < h; i++) {
1441
0
            lines = datas + i * wpls;
1442
0
            linef = dataf + i * wplf;
1443
0
            for (j = 2; j < w - 2; j++) {
1444
0
                val = GET_DATA_BYTE(lines, j - 2) +
1445
0
                      GET_DATA_BYTE(lines, j - 1) +
1446
0
                      GET_DATA_BYTE(lines, j) +
1447
0
                      GET_DATA_BYTE(lines, j + 1) +
1448
0
                      GET_DATA_BYTE(lines, j + 2);
1449
0
                linef[j] = val;
1450
0
            }
1451
0
        }
1452
0
    }
1453
1454
        /* Do vertical smoothing to finish the low-pass filter.
1455
         * At each pixel, if L is the lowpass value, I is the
1456
         * src pixel value and f is the fraction of highpass to
1457
         * be added to I, then the highpass filter value is
1458
         *     H = I - L
1459
         * and the new sharpened value is
1460
         *     N = I + f * H.                 */
1461
0
    if (halfwidth == 1) {
1462
0
        for (i = 1; i < h - 1; i++) {
1463
0
            linef0 = dataf + (i - 1) * wplf;
1464
0
            linef1 = dataf + i * wplf;
1465
0
            linef2 = dataf + (i + 1) * wplf;
1466
0
            lined = datad + i * wpld;
1467
0
            lines = datas + i * wpls;
1468
0
            norm = 1.0f / 9.0f;
1469
0
            for (j = 1; j < w - 1; j++) {
1470
0
                val = norm * (linef0[j] + linef1[j] +
1471
0
                              linef2[j]);         /* L: lowpass filter value */
1472
0
                sval = GET_DATA_BYTE(lines, j);   /* I: source pixel */
1473
0
                ival = (l_int32)(sval + fract * (sval - val) + 0.5);
1474
0
                ival = L_MAX(0, ival);
1475
0
                ival = L_MIN(255, ival);
1476
0
                SET_DATA_BYTE(lined, j, ival);
1477
0
            }
1478
0
        }
1479
0
    } else {
1480
0
        for (i = 2; i < h - 2; i++) {
1481
0
            linef0 = dataf + (i - 2) * wplf;
1482
0
            linef1 = dataf + (i - 1) * wplf;
1483
0
            linef2 = dataf + i * wplf;
1484
0
            linef3 = dataf + (i + 1) * wplf;
1485
0
            linef4 = dataf + (i + 2) * wplf;
1486
0
            lined = datad + i * wpld;
1487
0
            lines = datas + i * wpls;
1488
0
            norm = 1.0f / 25.0f;
1489
0
            for (j = 2; j < w - 2; j++) {
1490
0
                val = norm * (linef0[j] + linef1[j] + linef2[j] + linef3[j] +
1491
0
                              linef4[j]);  /* L: lowpass filter value */
1492
0
                sval = GET_DATA_BYTE(lines, j);   /* I: source pixel */
1493
0
                ival = (l_int32)(sval + fract * (sval - val) + 0.5);
1494
0
                ival = L_MAX(0, ival);
1495
0
                ival = L_MIN(255, ival);
1496
0
                SET_DATA_BYTE(lined, j, ival);
1497
0
            }
1498
0
        }
1499
0
    }
1500
1501
0
    fpixDestroy(&fpix);
1502
0
    return pixd;
1503
0
}
1504
1505
1506
/*-----------------------------------------------------------------------*
1507
 *                    Hue and saturation modification                    *
1508
 *-----------------------------------------------------------------------*/
1509
/*!
1510
 * \brief   pixModifyHue()
1511
 *
1512
 * \param[in]    pixd      [optional] can be null or equal to pixs
1513
 * \param[in]    pixs      32 bpp rgb
1514
 * \param[in]    fract     between -1.0 and 1.0
1515
 * \return  pixd, or NULL on error
1516
 *
1517
 * <pre>
1518
 * Notes:
1519
 *      (1) pixd must either be null or equal to pixs.
1520
 *          For in-place operation, set pixd == pixs:
1521
 *             pixEqualizeTRC(pixs, pixs, ...);
1522
 *          To get a new image, set pixd == null:
1523
 *             pixd = pixEqualizeTRC(NULL, pixs, ...);
1524
 *      (2) Use fract > 0.0 to increase hue value; < 0.0 to decrease it.
1525
 *          1.0 (or -1.0) represents a 360 degree rotation; i.e., no change.
1526
 *      (3) If no modification is requested (fract = -1.0 or 0 or 1.0),
1527
 *          return a copy unless in-place, in which case this is a no-op.
1528
 *      (4) This leaves saturation and intensity invariant.
1529
 *      (5) See discussion of color-modification methods, in coloring.c.
1530
 * </pre>
1531
 */
1532
PIX  *
1533
pixModifyHue(PIX       *pixd,
1534
             PIX       *pixs,
1535
             l_float32  fract)
1536
0
{
1537
0
l_int32    w, h, d, i, j, wpl, delhue;
1538
0
l_int32    rval, gval, bval, hval, sval, vval;
1539
0
l_uint32  *data, *line;
1540
1541
0
    if (!pixs)
1542
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1543
0
    if (pixGetColormap(pixs) != NULL)
1544
0
        return (PIX *)ERROR_PTR("pixs colormapped", __func__, NULL);
1545
0
    if (pixd && (pixd != pixs))
1546
0
        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
1547
0
    pixGetDimensions(pixs, &w, &h, &d);
1548
0
    if (d != 32)
1549
0
        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
1550
0
    if (L_ABS(fract) > 1.0)
1551
0
        return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", __func__, NULL);
1552
1553
0
    pixd = pixCopy(pixd, pixs);
1554
1555
0
    delhue = (l_int32)(240 * fract);
1556
0
    if (delhue == 0 || delhue == 240 || delhue == -240) {
1557
0
        L_WARNING("no change requested in hue\n", __func__);
1558
0
        return pixd;
1559
0
    }
1560
0
    if (delhue < 0)
1561
0
        delhue += 240;
1562
1563
0
    data = pixGetData(pixd);
1564
0
    wpl = pixGetWpl(pixd);
1565
0
    for (i = 0; i < h; i++) {
1566
0
        line = data + i * wpl;
1567
0
        for (j = 0; j < w; j++) {
1568
0
            extractRGBValues(line[j], &rval, &gval, &bval);
1569
0
            convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
1570
0
            hval = (hval + delhue) % 240;
1571
0
            convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
1572
0
            composeRGBPixel(rval, gval, bval, line + j);
1573
0
        }
1574
0
    }
1575
0
    if (pixGetSpp(pixs) == 4)
1576
0
        pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
1577
1578
0
    return pixd;
1579
0
}
1580
1581
1582
/*!
1583
 * \brief   pixModifySaturation()
1584
 *
1585
 * \param[in]    pixd     [optional] can be null, existing or equal to pixs
1586
 * \param[in]    pixs     32 bpp rgb
1587
 * \param[in]    fract    between -1.0 and 1.0
1588
 * \return  pixd, or NULL on error
1589
 *
1590
 * <pre>
1591
 * Notes:
1592
 *      (1) If fract > 0.0, it gives the fraction that the pixel
1593
 *          saturation is moved from its initial value toward 255.
1594
 *          If fract < 0.0, it gives the fraction that the pixel
1595
 *          saturation is moved from its initial value toward 0.
1596
 *          The limiting values for fract = -1.0 (1.0) thus set the
1597
 *          saturation to 0 (255).
1598
 *      (2) If fract = 0, no modification is requested; return a copy
1599
 *          unless in-place, in which case this is a no-op.
1600
 *      (3) This leaves hue and intensity invariant.
1601
 *      (4) See discussion of color-modification methods, in coloring.c.
1602
 * </pre>
1603
 */
1604
PIX  *
1605
pixModifySaturation(PIX       *pixd,
1606
                    PIX       *pixs,
1607
                    l_float32  fract)
1608
0
{
1609
0
l_int32    w, h, d, i, j, wpl;
1610
0
l_int32    rval, gval, bval, hval, sval, vval;
1611
0
l_uint32  *data, *line;
1612
1613
0
    if (!pixs)
1614
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1615
0
    pixGetDimensions(pixs, &w, &h, &d);
1616
0
    if (d != 32)
1617
0
        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
1618
0
    if (L_ABS(fract) > 1.0)
1619
0
        return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", __func__, NULL);
1620
1621
0
    pixd = pixCopy(pixd, pixs);
1622
0
    if (fract == 0.0) {
1623
0
        L_WARNING("no change requested in saturation\n", __func__);
1624
0
        return pixd;
1625
0
    }
1626
1627
0
    data = pixGetData(pixd);
1628
0
    wpl = pixGetWpl(pixd);
1629
0
    for (i = 0; i < h; i++) {
1630
0
        line = data + i * wpl;
1631
0
        for (j = 0; j < w; j++) {
1632
0
            extractRGBValues(line[j], &rval, &gval, &bval);
1633
0
            convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
1634
0
            if (fract < 0.0)
1635
0
                sval = (l_int32)(sval * (1.0 + fract));
1636
0
            else
1637
0
                sval = (l_int32)(sval + fract * (255 - sval));
1638
0
            convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
1639
0
            composeRGBPixel(rval, gval, bval, line + j);
1640
0
        }
1641
0
    }
1642
0
    if (pixGetSpp(pixs) == 4)
1643
0
        pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
1644
1645
0
    return pixd;
1646
0
}
1647
1648
1649
/*!
1650
 * \brief   pixMeasureSaturation()
1651
 *
1652
 * \param[in]    pixs     32 bpp rgb
1653
 * \param[in]    factor   subsampling factor; integer >= 1
1654
 * \param[out]   psat     average saturation
1655
 * \return  0 if OK, 1 on error
1656
 */
1657
l_int32
1658
pixMeasureSaturation(PIX        *pixs,
1659
                     l_int32     factor,
1660
                     l_float32  *psat)
1661
0
{
1662
0
l_int32    w, h, d, i, j, wpl, sum, count;
1663
0
l_int32    rval, gval, bval, hval, sval, vval;
1664
0
l_uint32  *data, *line;
1665
1666
0
    if (!psat)
1667
0
        return ERROR_INT("pixs not defined", __func__, 1);
1668
0
    *psat = 0.0;
1669
0
    if (!pixs)
1670
0
        return ERROR_INT("pixs not defined", __func__, 1);
1671
0
    pixGetDimensions(pixs, &w, &h, &d);
1672
0
    if (d != 32)
1673
0
        return ERROR_INT("pixs not 32 bpp", __func__, 1);
1674
0
    if (factor < 1)
1675
0
        return ERROR_INT("subsampling factor < 1", __func__, 1);
1676
1677
0
    data = pixGetData(pixs);
1678
0
    wpl = pixGetWpl(pixs);
1679
0
    for (i = 0, sum = 0, count = 0; i < h; i += factor) {
1680
0
        line = data + i * wpl;
1681
0
        for (j = 0; j < w; j += factor) {
1682
0
            extractRGBValues(line[j], &rval, &gval, &bval);
1683
0
            convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
1684
0
            sum += sval;
1685
0
            count++;
1686
0
        }
1687
0
    }
1688
1689
0
    if (count > 0)
1690
0
        *psat = (l_float32)sum / (l_float32)count;
1691
0
    return 0;
1692
0
}
1693
1694
1695
/*!
1696
 * \brief   pixModifyBrightness()
1697
 *
1698
 * \param[in]    pixd     [optional] can be null, existing or equal to pixs
1699
 * \param[in]    pixs     32 bpp rgb
1700
 * \param[in]    fract    between -1.0 and 1.0
1701
 * \return  pixd, or NULL on error
1702
 *
1703
 * <pre>
1704
 * Notes:
1705
 *      (1) If fract > 0.0, it gives the fraction that the v-parameter,
1706
 *          which is max(r,g,b), is moved from its initial value toward 255.
1707
 *          If fract < 0.0, it gives the fraction that the v-parameter
1708
 *          is moved from its initial value toward 0.
1709
 *          The limiting values for fract = -1.0 (1.0) thus set the
1710
 *          v-parameter to 0 (255).
1711
 *      (2) If fract = 0, no modification is requested; return a copy
1712
 *          unless in-place, in which case this is a no-op.
1713
 *      (3) This leaves hue and saturation invariant.
1714
 *      (4) See discussion of color-modification methods, in coloring.c.
1715
 * </pre>
1716
 */
1717
PIX  *
1718
pixModifyBrightness(PIX       *pixd,
1719
                    PIX       *pixs,
1720
                    l_float32  fract)
1721
0
{
1722
0
l_int32    w, h, d, i, j, wpl;
1723
0
l_int32    rval, gval, bval, hval, sval, vval;
1724
0
l_uint32  *data, *line;
1725
1726
0
    if (!pixs)
1727
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1728
0
    pixGetDimensions(pixs, &w, &h, &d);
1729
0
    if (d != 32)
1730
0
        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
1731
0
    if (L_ABS(fract) > 1.0)
1732
0
        return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", __func__, NULL);
1733
1734
0
    pixd = pixCopy(pixd, pixs);
1735
0
    if (fract == 0.0) {
1736
0
        L_WARNING("no change requested in brightness\n", __func__);
1737
0
        return pixd;
1738
0
    }
1739
1740
0
    data = pixGetData(pixd);
1741
0
    wpl = pixGetWpl(pixd);
1742
0
    for (i = 0; i < h; i++) {
1743
0
        line = data + i * wpl;
1744
0
        for (j = 0; j < w; j++) {
1745
0
            extractRGBValues(line[j], &rval, &gval, &bval);
1746
0
            convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
1747
0
            if (fract > 0.0)
1748
0
                vval = (l_int32)(vval + fract * (255.0 - vval));
1749
0
            else
1750
0
                vval = (l_int32)(vval * (1.0 + fract));
1751
0
            convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
1752
0
            composeRGBPixel(rval, gval, bval, line + j);
1753
0
        }
1754
0
    }
1755
0
    if (pixGetSpp(pixs) == 4)
1756
0
        pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
1757
1758
0
    return pixd;
1759
0
}
1760
1761
1762
/*-----------------------------------------------------------------------*
1763
 *                             Color shifting                            *
1764
 *-----------------------------------------------------------------------*/
1765
/*!
1766
 * \brief   pixMosaicColorShiftRGB()
1767
 *
1768
 * \param[in]    pixs     32 bpp rgb
1769
 * \param[in]    roff   center offset of red component
1770
 * \param[in]    goff   center offset of green component
1771
 * \param[in]    boff   center offset of blue component
1772
 * \param[in]    delta  increments from center offsets [0.0 - 0.1];
1773
 *                      use 0.0 to get the default (0.04)
1774
 * \param[in]    nincr  number of increments in each (positive and negative)
1775
 *                      direction; use 0 to get the default (2).
1776
 * \return  pix, or NULL on error
1777
 *
1778
 * <pre>
1779
 * Notes:
1780
 *      (1) This generates a mosaic view of the effect of shifting the RGB
1781
 *          components.  See pixColorShiftRGB() for details on the shifting.
1782
 *      (2) The offsets (%roff, %goff, %boff) set the color center point,
1783
 *          and the deviations from this are shown separately for deltas
1784
 *          in r, g and b.  For each component, we show 2 * %nincr + 1
1785
 *          images.
1786
 *      (3) The pix must have minimum dimensions of 100 and an aspect
1787
 *          ratio not exceeding 5.0.
1788
 *      (4) Usage: color prints differ from the original due to three factors:
1789
 *          illumination, calibration of the camera in acquisition,
1790
 *          and calibration of the printer.  This function can be used
1791
 *          to iteratively match a color print to the original.  On each
1792
 *          iteration, the center offsets are set to the best match so
1793
 *          far, and the %delta increments are typically reduced.
1794
 * </pre>
1795
 */
1796
PIX *
1797
pixMosaicColorShiftRGB(PIX       *pixs,
1798
                       l_float32  roff,
1799
                       l_float32  goff,
1800
                       l_float32  boff,
1801
                       l_float32  delta,
1802
                       l_int32    nincr)
1803
0
{
1804
0
char       buf[64];
1805
0
l_int32    i, w, h;
1806
0
l_float32  del, ratio;
1807
0
L_BMF     *bmf;
1808
0
PIX       *pix1, *pix2, *pix3;
1809
0
PIXA      *pixa;
1810
1811
0
    if (!pixs  || pixGetDepth(pixs) != 32)
1812
0
        return (PIX *)ERROR_PTR("pixs undefined or not rgb", __func__, NULL);
1813
0
    if (roff < -1.0 || roff > 1.0)
1814
0
        return (PIX *)ERROR_PTR("roff not in [-1.0, 1.0]", __func__, NULL);
1815
0
    if (goff < -1.0 || goff > 1.0)
1816
0
        return (PIX *)ERROR_PTR("goff not in [-1.0, 1.0]", __func__, NULL);
1817
0
    if (boff < -1.0 || boff > 1.0)
1818
0
        return (PIX *)ERROR_PTR("boff not in [-1.0, 1.0]", __func__, NULL);
1819
0
    if (delta < 0.0 || delta > 0.1)
1820
0
        return (PIX *)ERROR_PTR("delta not in [0.0, 0.1]", __func__, NULL);
1821
0
    if (delta == 0.0) delta = 0.04f;
1822
0
    if (nincr < 0 || nincr > 6)
1823
0
        return (PIX *)ERROR_PTR("nincr not in [0, 6]", __func__, NULL);
1824
0
    if (nincr == 0) nincr = 2;
1825
1826
        /* Require width and height to be >= 100, and the aspect ratio <= 5.0 */
1827
0
    pixGetDimensions(pixs, &w, &h, NULL);
1828
0
    if (w < 100 || h < 100)
1829
0
        return (PIX *)ERROR_PTR("w and h not both >= 100", __func__, NULL);
1830
0
    pixMaxAspectRatio(pixs, &ratio);
1831
0
    if (ratio < 1.0 || ratio > 5.0) {
1832
0
        L_ERROR("invalid aspect ratio %5.1f\n", __func__, ratio);
1833
0
        return NULL;
1834
0
    }
1835
1836
0
    pixa = pixaCreate(3 * (2 * nincr + 1));
1837
0
    bmf = bmfCreate(NULL, 8);
1838
0
    pix1 = pixScaleToSize(pixs, 400, 0);
1839
0
    for (i = 0, del = - nincr * delta; i < 2 * nincr + 1; i++, del += delta) {
1840
0
        pix2 = pixColorShiftRGB(pix1, roff + del, goff, boff);
1841
0
        snprintf(buf, sizeof(buf), "%4.2f, %4.2f, %4.2f",
1842
0
                 roff + del, goff, boff);
1843
0
        pix3 = pixAddSingleTextblock(pix2, bmf, buf, 0xff000000,
1844
0
                                     L_ADD_BELOW, 0);
1845
0
        pixaAddPix(pixa, pix3, L_INSERT);
1846
0
        pixDestroy(&pix2);
1847
0
    }
1848
0
    for (i = 0, del = - nincr * delta; i < 2 * nincr + 1; i++, del += delta) {
1849
0
        pix2 = pixColorShiftRGB(pix1, roff, goff + del, boff);
1850
0
        snprintf(buf, sizeof(buf), "%4.2f, %4.2f, %4.2f",
1851
0
                 roff, goff + del, boff);
1852
0
        pix3 = pixAddSingleTextblock(pix2, bmf, buf, 0xff000000,
1853
0
                                     L_ADD_BELOW, 0);
1854
0
        pixaAddPix(pixa, pix3, L_INSERT);
1855
0
        pixDestroy(&pix2);
1856
0
    }
1857
0
    for (i = 0, del = - nincr * delta; i < 2 * nincr + 1; i++, del += delta) {
1858
0
        pix2 = pixColorShiftRGB(pix1, roff, goff, boff + del);
1859
0
        snprintf(buf, sizeof(buf), "%4.2f, %4.2f, %4.2f",
1860
0
                 roff, goff, boff + del);
1861
0
        pix3 = pixAddSingleTextblock(pix2, bmf, buf, 0xff000000,
1862
0
                                     L_ADD_BELOW, 0);
1863
0
        pixaAddPix(pixa, pix3, L_INSERT);
1864
0
        pixDestroy(&pix2);
1865
0
    }
1866
0
    pixDestroy(&pix1);
1867
1868
0
    pix1 = pixaDisplayTiledAndScaled(pixa, 32, 300, 2 * nincr + 1, 0, 30, 2);
1869
0
    pixaDestroy(&pixa);
1870
0
    bmfDestroy(&bmf);
1871
0
    return pix1;
1872
0
}
1873
1874
1875
/*!
1876
 * \brief   pixColorShiftRGB()
1877
 *
1878
 * \param[in]    pixs     32 bpp rgb
1879
 * \param[in]    rfract   fractional shift in red component
1880
 * \param[in]    gfract   fractional shift in green component
1881
 * \param[in]    bfract   fractional shift in blue component
1882
 * \return  pixd, or NULL on error
1883
 *
1884
 * <pre>
1885
 * Notes:
1886
 *      (1) This allows independent fractional shifts of the r,g and b
1887
 *          components.  A positive shift pushes to saturation (255);
1888
 *          a negative shift pushes toward 0 (black).
1889
 *      (2) The effect can be imagined using a color wheel that consists
1890
 *          (for our purposes) of these 6 colors, separated by 60 degrees:
1891
 *             red, magenta, blue, cyan, green, yellow
1892
 *      (3) So, for example, a negative shift of the blue component
1893
 *          (bfract < 0) could be accompanied by positive shifts
1894
 *          of red and green to make an image more yellow.
1895
 *      (4) Examples of limiting cases:
1896
 *            rfract = 1 ==> r = 255
1897
 *            rfract = -1 ==> r = 0
1898
 * </pre>
1899
 */
1900
PIX *
1901
pixColorShiftRGB(PIX       *pixs,
1902
                 l_float32  rfract,
1903
                 l_float32  gfract,
1904
                 l_float32  bfract)
1905
0
{
1906
0
l_int32    w, h, i, j, wpls, wpld, rval, gval, bval;
1907
0
l_int32   *rlut, *glut, *blut;
1908
0
l_uint32  *datas, *datad, *lines, *lined;
1909
0
l_float32  fi;
1910
0
PIX       *pixd;
1911
1912
0
    if (!pixs)
1913
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1914
0
    if (pixGetDepth(pixs) != 32)
1915
0
        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
1916
0
    if (rfract < -1.0 || rfract > 1.0)
1917
0
        return (PIX *)ERROR_PTR("rfract not in [-1.0, 1.0]", __func__, NULL);
1918
0
    if (gfract < -1.0 || gfract > 1.0)
1919
0
        return (PIX *)ERROR_PTR("gfract not in [-1.0, 1.0]", __func__, NULL);
1920
0
    if (bfract < -1.0 || bfract > 1.0)
1921
0
        return (PIX *)ERROR_PTR("bfract not in [-1.0, 1.0]", __func__, NULL);
1922
0
    if (rfract == 0.0 && gfract == 0.0 && bfract == 0.0)
1923
0
        return pixCopy(NULL, pixs);
1924
1925
0
    rlut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
1926
0
    glut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
1927
0
    blut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
1928
0
    for (i = 0; i < 256; i++) {
1929
0
        fi = i;
1930
0
        if (rfract >= 0) {
1931
0
            rlut[i] = (l_int32)(fi + (255.0 - fi) * rfract);
1932
0
        } else {
1933
0
            rlut[i] = (l_int32)(fi * (1.0 + rfract));
1934
0
        }
1935
0
        if (gfract >= 0) {
1936
0
            glut[i] = (l_int32)(fi + (255.0 - fi) * gfract);
1937
0
        } else {
1938
0
            glut[i] = (l_int32)(fi * (1.0 + gfract));
1939
0
        }
1940
0
        if (bfract >= 0) {
1941
0
            blut[i] = (l_int32)(fi + (255.0 - fi) * bfract);
1942
0
        } else {
1943
0
            blut[i] = (l_int32)(fi * (1.0 + bfract));
1944
0
        }
1945
0
    }
1946
1947
0
    pixGetDimensions(pixs, &w, &h, NULL);
1948
0
    datas = pixGetData(pixs);
1949
0
    wpls = pixGetWpl(pixs);
1950
0
    pixd = pixCreate(w, h, 32);
1951
0
    datad = pixGetData(pixd);
1952
0
    wpld = pixGetWpl(pixd);
1953
0
    for (i = 0; i < h; i++) {
1954
0
        lines = datas + i * wpls;
1955
0
        lined = datad + i * wpld;
1956
0
        for (j = 0; j < w; j++) {
1957
0
            extractRGBValues(lines[j], &rval, &gval, &bval);
1958
0
            composeRGBPixel(rlut[rval], glut[gval], blut[bval], lined + j);
1959
0
        }
1960
0
    }
1961
1962
0
    LEPT_FREE(rlut);
1963
0
    LEPT_FREE(glut);
1964
0
    LEPT_FREE(blut);
1965
0
    return pixd;
1966
0
}
1967
1968
/*-----------------------------------------------------------------------*
1969
 *                     Darken gray (unsaturated) pixels
1970
 *-----------------------------------------------------------------------*/
1971
/*!
1972
 * \brief   pixDarkenGray()
1973
 *
1974
 * \param[in]    pixd      [optional] can be null or equal to pixs
1975
 * \param[in]    pixs      32 bpp rgb
1976
 * \param[in]    thresh    pixels with max component >= %thresh are unchanged
1977
 * \param[in]    satlimit  pixels with saturation >= %satlimit are unchanged
1978
 * \return  pixd, or NULL on error
1979
 *
1980
 * <pre>
1981
 * Notes:
1982
 *      (1) This darkens gray pixels, by a fraction (sat/%satlimit), where
1983
 *          the saturation, sat, is the component difference (max - min).
1984
 *          The pixel value is unchanged if sat >= %satlimit.  A typical
1985
 *          value of %satlimit might be 40; the larger the value, the
1986
 *          more that pixels with a smaller saturation will be darkened.
1987
 *      (2) Pixels with max component >= %thresh are unchanged. This can be
1988
 *          used to prevent bright pixels with low saturation from being
1989
 *          darkened.  Setting thresh == 0 is a no-op; setting %thresh == 255
1990
 *          causes the darkening to be applied to all pixels.
1991
 *      (3) This function is useful to enhance pixels relative to a
1992
 *          gray background.
1993
 *      (4) A related function that builds a 1 bpp mask over the gray
1994
 *          pixels is pixMaskOverGrayPixels().
1995
 * </pre>
1996
 */
1997
PIX *
1998
pixDarkenGray(PIX     *pixd,
1999
              PIX     *pixs,
2000
              l_int32  thresh,
2001
              l_int32  satlimit)
2002
0
{
2003
0
l_int32    w, h, i, j, wpls, wpld;
2004
0
l_int32    rval, gval, bval, minrg, min, maxrg, max, sat;
2005
0
l_uint32  *datas, *datad, *lines, *lined;
2006
0
l_float32  ratio;
2007
2008
0
    if (!pixs || pixGetDepth(pixs) != 32)
2009
0
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", __func__, NULL);
2010
0
    if (thresh < 0 || thresh > 255)
2011
0
        return (PIX *)ERROR_PTR("invalid thresh", __func__, NULL);
2012
0
    if (satlimit < 1)
2013
0
        return (PIX *)ERROR_PTR("invalid satlimit", __func__, NULL);
2014
0
    if (pixd && (pixs != pixd))
2015
0
        return (PIX *)ERROR_PTR("not new or in-place", __func__, NULL);
2016
2017
0
    pixGetDimensions(pixs, &w, &h, NULL);
2018
0
    datas = pixGetData(pixs);
2019
0
    wpls = pixGetWpl(pixs);
2020
0
    if ((pixd = pixCopy(pixd, pixs)) == NULL)
2021
0
        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
2022
0
    datad = pixGetData(pixd);
2023
0
    wpld = pixGetWpl(pixd);
2024
2025
0
    for (i = 0; i < h; i++) {
2026
0
        lines = datas + i * wpls;
2027
0
        lined = datad + i * wpld;
2028
0
        for (j = 0; j < w; j++) {
2029
0
            extractRGBValues(lines[j], &rval, &gval, &bval);
2030
0
            minrg = L_MIN(rval, gval);
2031
0
            min = L_MIN(minrg, bval);
2032
0
            maxrg = L_MAX(rval, gval);
2033
0
            max = L_MAX(maxrg, bval);
2034
0
            sat = max - min;
2035
0
            if (max >= thresh || sat >= satlimit)
2036
0
                continue;
2037
0
            ratio = (l_float32)sat / (l_float32)satlimit;
2038
0
            composeRGBPixel((l_int32)(ratio * rval), (l_int32)(ratio * gval),
2039
0
                            (l_int32)(ratio * bval), &lined[j]);
2040
0
        }
2041
0
    }
2042
0
    return pixd;
2043
0
}
2044
2045
2046
/*-----------------------------------------------------------------------*
2047
 *            General multiplicative constant color transform            *
2048
 *-----------------------------------------------------------------------*/
2049
/*!
2050
 * \brief   pixMultConstantColor()
2051
 *
2052
 * \param[in]    pixs     colormapped or rgb
2053
 * \param[in]    rfact    red multiplicative factor
2054
 * \param[in]    gfact    green multiplicative factor
2055
 * \param[in]    bfact    blue multiplicative factor
2056
 * \return  pixd colormapped or rgb, with colors scaled, or NULL on error
2057
 *
2058
 * <pre>
2059
 * Notes:
2060
 *      (1) rfact, gfact and bfact can only have non-negative values.
2061
 *          They can be greater than 1.0.  All transformed component
2062
 *          values are clipped to the interval [0, 255].
2063
 *      (2) For multiplication with a general 3x3 matrix of constants,
2064
 *          use pixMultMatrixColor().
2065
 * </pre>
2066
 */
2067
PIX *
2068
pixMultConstantColor(PIX       *pixs,
2069
                     l_float32  rfact,
2070
                     l_float32  gfact,
2071
                     l_float32  bfact)
2072
0
{
2073
0
l_int32    i, j, w, h, d, wpls, wpld;
2074
0
l_int32    ncolors, rval, gval, bval, nrval, ngval, nbval;
2075
0
l_uint32   nval;
2076
0
l_uint32  *datas, *datad, *lines, *lined;
2077
0
PIX       *pixd;
2078
0
PIXCMAP   *cmap;
2079
2080
0
    if (!pixs)
2081
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
2082
0
    pixGetDimensions(pixs, &w, &h, &d);
2083
0
    cmap = pixGetColormap(pixs);
2084
0
    if (!cmap && d != 32)
2085
0
        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", __func__, NULL);
2086
0
    rfact = L_MAX(0.0, rfact);
2087
0
    gfact = L_MAX(0.0, gfact);
2088
0
    bfact = L_MAX(0.0, bfact);
2089
2090
0
    if (cmap) {
2091
0
        if ((pixd = pixCopy(NULL, pixs)) == NULL)
2092
0
            return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
2093
0
        cmap = pixGetColormap(pixd);
2094
0
        ncolors = pixcmapGetCount(cmap);
2095
0
        for (i = 0; i < ncolors; i++) {
2096
0
            pixcmapGetColor(cmap, i, &rval, &gval, &bval);
2097
0
            nrval = (l_int32)(rfact * rval);
2098
0
            ngval = (l_int32)(gfact * gval);
2099
0
            nbval = (l_int32)(bfact * bval);
2100
0
            nrval = L_MIN(255, nrval);
2101
0
            ngval = L_MIN(255, ngval);
2102
0
            nbval = L_MIN(255, nbval);
2103
0
            pixcmapResetColor(cmap, i, nrval, ngval, nbval);
2104
0
        }
2105
0
        return pixd;
2106
0
    }
2107
2108
0
    if ((pixd = pixCreateTemplate(pixs)) == NULL)
2109
0
        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
2110
0
    datas = pixGetData(pixs);
2111
0
    datad = pixGetData(pixd);
2112
0
    wpls = pixGetWpl(pixs);
2113
0
    wpld = pixGetWpl(pixd);
2114
0
    for (i = 0; i < h; i++) {
2115
0
        lines = datas + i * wpls;
2116
0
        lined = datad + i * wpld;
2117
0
        for (j = 0; j < w; j++) {
2118
0
            extractRGBValues(lines[j], &rval, &gval, &bval);
2119
0
            nrval = (l_int32)(rfact * rval);
2120
0
            ngval = (l_int32)(gfact * gval);
2121
0
            nbval = (l_int32)(bfact * bval);
2122
0
            nrval = L_MIN(255, nrval);
2123
0
            ngval = L_MIN(255, ngval);
2124
0
            nbval = L_MIN(255, nbval);
2125
0
            composeRGBPixel(nrval, ngval, nbval, &nval);
2126
0
            *(lined + j) = nval;
2127
0
        }
2128
0
    }
2129
2130
0
    return pixd;
2131
0
}
2132
2133
2134
/*!
2135
 * \brief   pixMultMatrixColor()
2136
 *
2137
 * \param[in]    pixs    colormapped or rgb
2138
 * \param[in]    kel     kernel 3x3 matrix of floats
2139
 * \return  pixd colormapped or rgb, or NULL on error
2140
 *
2141
 * <pre>
2142
 * Notes:
2143
 *      (1) The kernel is a data structure used mostly for floating point
2144
 *          convolution.  Here it is a 3x3 matrix of floats that are used
2145
 *          to transform the pixel values by matrix multiplication:
2146
 *            nrval = a[0,0] * rval + a[0,1] * gval + a[0,2] * bval
2147
 *            ngval = a[1,0] * rval + a[1,1] * gval + a[1,2] * bval
2148
 *            nbval = a[2,0] * rval + a[2,1] * gval + a[2,2] * bval
2149
 *      (2) The matrix can be generated in several ways.
2150
 *          See kernel.c for details.  Here are two of them:
2151
 *            (a) kel = kernelCreate(3, 3);
2152
 *                kernelSetElement(kel, 0, 0, val00);
2153
 *                kernelSetElement(kel, 0, 1, val01);
2154
 *                ...
2155
 *            (b) from a static string; e.g.,:
2156
 *                const char *kdata = " 0.6  0.3 -0.2 "
2157
 *                                    " 0.1  1.2  0.4 "
2158
 *                                    " -0.4 0.2  0.9 ";
2159
 *                kel = kernelCreateFromString(3, 3, 0, 0, kdata);
2160
 *      (3) For the special case where the matrix is diagonal, it is easier
2161
 *          to use pixMultConstantColor().
2162
 *      (4) Matrix entries can have positive and negative values, and can
2163
 *          be larger than 1.0.  All transformed component values
2164
 *          are clipped to [0, 255].
2165
 * </pre>
2166
 */
2167
PIX *
2168
pixMultMatrixColor(PIX       *pixs,
2169
                   L_KERNEL  *kel)
2170
0
{
2171
0
l_int32    i, j, index, kw, kh, w, h, d, wpls, wpld;
2172
0
l_int32    ncolors, rval, gval, bval, nrval, ngval, nbval;
2173
0
l_uint32   nval;
2174
0
l_uint32  *datas, *datad, *lines, *lined;
2175
0
l_float32  v[9];  /* use linear array for convenience */
2176
0
PIX       *pixd;
2177
0
PIXCMAP   *cmap;
2178
2179
0
    if (!pixs)
2180
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
2181
0
    if (!kel)
2182
0
        return (PIX *)ERROR_PTR("kel not defined", __func__, NULL);
2183
0
    kernelGetParameters(kel, &kw, &kh, NULL, NULL);
2184
0
    if (kw != 3 || kh != 3)
2185
0
        return (PIX *)ERROR_PTR("matrix not 3x3", __func__, NULL);
2186
0
    pixGetDimensions(pixs, &w, &h, &d);
2187
0
    cmap = pixGetColormap(pixs);
2188
0
    if (!cmap && d != 32)
2189
0
        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", __func__, NULL);
2190
2191
0
    for (i = 0, index = 0; i < 3; i++)
2192
0
        for (j = 0; j < 3; j++, index++)
2193
0
            kernelGetElement(kel, i, j, v + index);
2194
2195
0
    if (cmap) {
2196
0
        if ((pixd = pixCopy(NULL, pixs)) == NULL)
2197
0
            return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
2198
0
        cmap = pixGetColormap(pixd);
2199
0
        ncolors = pixcmapGetCount(cmap);
2200
0
        for (i = 0; i < ncolors; i++) {
2201
0
            pixcmapGetColor(cmap, i, &rval, &gval, &bval);
2202
0
            nrval = (l_int32)(v[0] * rval + v[1] * gval + v[2] * bval);
2203
0
            ngval = (l_int32)(v[3] * rval + v[4] * gval + v[5] * bval);
2204
0
            nbval = (l_int32)(v[6] * rval + v[7] * gval + v[8] * bval);
2205
0
            nrval = L_MAX(0, L_MIN(255, nrval));
2206
0
            ngval = L_MAX(0, L_MIN(255, ngval));
2207
0
            nbval = L_MAX(0, L_MIN(255, nbval));
2208
0
            pixcmapResetColor(cmap, i, nrval, ngval, nbval);
2209
0
        }
2210
0
        return pixd;
2211
0
    }
2212
2213
0
    if ((pixd = pixCreateTemplate(pixs)) == NULL)
2214
0
        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
2215
0
    datas = pixGetData(pixs);
2216
0
    datad = pixGetData(pixd);
2217
0
    wpls = pixGetWpl(pixs);
2218
0
    wpld = pixGetWpl(pixd);
2219
0
    for (i = 0; i < h; i++) {
2220
0
        lines = datas + i * wpls;
2221
0
        lined = datad + i * wpld;
2222
0
        for (j = 0; j < w; j++) {
2223
0
            extractRGBValues(lines[j], &rval, &gval, &bval);
2224
0
            nrval = (l_int32)(v[0] * rval + v[1] * gval + v[2] * bval);
2225
0
            ngval = (l_int32)(v[3] * rval + v[4] * gval + v[5] * bval);
2226
0
            nbval = (l_int32)(v[6] * rval + v[7] * gval + v[8] * bval);
2227
0
            nrval = L_MAX(0, L_MIN(255, nrval));
2228
0
            ngval = L_MAX(0, L_MIN(255, ngval));
2229
0
            nbval = L_MAX(0, L_MIN(255, nbval));
2230
0
            composeRGBPixel(nrval, ngval, nbval, &nval);
2231
0
            *(lined + j) = nval;
2232
0
        }
2233
0
    }
2234
2235
0
    return pixd;
2236
0
}
2237
2238
2239
/*-------------------------------------------------------------*
2240
 *                    Half-edge by bandpass                    *
2241
 *-------------------------------------------------------------*/
2242
/*!
2243
 * \brief   pixHalfEdgeByBandpass()
2244
 *
2245
 * \param[in]    pixs         8 bpp gray or 32 bpp rgb
2246
 * \param[in]    sm1h, sm1v   "half-widths" of smoothing filter sm1
2247
 * \param[in]    sm2h, sm2v   "half-widths" of smoothing filter sm2;
2248
 *                            require sm2 != sm1
2249
 * \return  pixd, or NULL on error
2250
 *
2251
 * <pre>
2252
 * Notes:
2253
 *      (1) We use symmetric smoothing filters of odd dimension,
2254
 *          typically use 3, 5, 7, etc.  The smoothing parameters
2255
 *          for these are 1, 2, 3, etc.  The filter size is related
2256
 *          to the smoothing parameter by
2257
 *               size = 2 * smoothing + 1
2258
 *      (2) Because we take the difference of two lowpass filters,
2259
 *          this is actually a bandpass filter.
2260
 *      (3) We allow both filters to be anisotropic.
2261
 *      (4) Consider either the h or v component of the 2 filters.
2262
 *          Depending on whether sm1 > sm2 or sm2 > sm1, we get
2263
 *          different halves of the smoothed gradients (or "edges").
2264
 *          This difference of smoothed signals looks more like
2265
 *          a second derivative of a transition, which we rectify
2266
 *          by not allowing the signal to go below zero.  If sm1 < sm2,
2267
 *          the sm2 transition is broader, so the difference between
2268
 *          sm1 and sm2 signals is positive on the upper half of
2269
 *          the transition.  Likewise, if sm1 > sm2, the sm1 - sm2
2270
 *          signal difference is positive on the lower half of
2271
 *          the transition.
2272
 * </pre>
2273
 */
2274
PIX *
2275
pixHalfEdgeByBandpass(PIX     *pixs,
2276
                      l_int32  sm1h,
2277
                      l_int32  sm1v,
2278
                      l_int32  sm2h,
2279
                      l_int32  sm2v)
2280
0
{
2281
0
l_int32  d;
2282
0
PIX     *pixg, *pixacc, *pixc1, *pixc2;
2283
2284
0
    if (!pixs)
2285
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
2286
0
    if (sm1h == sm2h && sm1v == sm2v)
2287
0
        return (PIX *)ERROR_PTR("sm2 = sm1", __func__, NULL);
2288
0
    d = pixGetDepth(pixs);
2289
0
    if (d != 8 && d != 32)
2290
0
        return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL);
2291
0
    if (d == 32)
2292
0
        pixg = pixConvertRGBToLuminance(pixs);
2293
0
    else   /* d == 8 */
2294
0
        pixg = pixClone(pixs);
2295
2296
        /* Make a convolution accumulator and use it twice */
2297
0
    if ((pixacc = pixBlockconvAccum(pixg)) == NULL) {
2298
0
        pixDestroy(&pixg);
2299
0
        return (PIX *)ERROR_PTR("pixacc not made", __func__, NULL);
2300
0
    }
2301
0
    if ((pixc1 = pixBlockconvGray(pixg, pixacc, sm1h, sm1v)) == NULL) {
2302
0
        pixDestroy(&pixg);
2303
0
        pixDestroy(&pixacc);
2304
0
        return (PIX *)ERROR_PTR("pixc1 not made", __func__, NULL);
2305
0
    }
2306
0
    pixc2 = pixBlockconvGray(pixg, pixacc, sm2h, sm2v);
2307
0
    pixDestroy(&pixg);
2308
0
    pixDestroy(&pixacc);
2309
0
    if (!pixc2) {
2310
0
        pixDestroy(&pixc1);
2311
0
        return (PIX *)ERROR_PTR("pixc2 not made", __func__, NULL);
2312
0
    }
2313
2314
        /* Compute the half-edge using pixc1 - pixc2.  */
2315
0
    pixSubtractGray(pixc1, pixc1, pixc2);
2316
0
    pixDestroy(&pixc2);
2317
0
    return pixc1;
2318
0
}