1 /*====================================================================*
2 - Copyright (C) 2001 Leptonica. All rights reserved.
3 - This software is distributed in the hope that it will be
4 - useful, but with NO WARRANTY OF ANY KIND.
5 - No author or distributor accepts responsibility to anyone for the
6 - consequences of using this software, or for whether it serves any
7 - particular purpose or works at all, unless he or she says so in
8 - writing. Everyone is granted permission to copy, modify and
9 - redistribute this source code, for commercial or non-commercial
10 - purposes, with the following restrictions: (1) the origin of this
11 - source code must not be misrepresented; (2) modified versions must
12 - be plainly marked as such; and (3) this notice may not be removed
13 - or altered from any source or modified source distribution.
14 *====================================================================*/
15
16
17 /*
18 * pixarith.c
19 *
20 * One-image grayscale arithmetic operations (8, 16, 32 bpp)
21 * l_int32 pixAddConstantGray()
22 * l_int32 pixMultConstantGray()
23 *
24 * Two-image grayscale arithmetic operations (8, 16, 32 bpp)
25 * PIX *pixAddGray()
26 * PIX *pixSubtractGray()
27 *
28 * Grayscale threshold operation (8, 16, 32 bpp)
29 * PIX *pixThresholdToValue()
30 *
31 * Image accumulator arithmetic operations
32 * PIX *pixInitAccumulate()
33 * PIX *pixFinalAccumulate()
34 * PIX *pixFinalAccumulateThreshold()
35 * l_int32 pixAccumulate()
36 * l_int32 pixMultConstAccumulate()
37 *
38 * Absolute value of difference
39 * PIX *pixAbsDifference()
40 *
41 * Two-image min and max operations (8 and 16 bpp)
42 * PIX *pixMinOrMax()
43 *
44 * Scale pix for maximum dynamic range in 8 bpp image:
45 * PIX *pixMaxDynamicRange()
46 *
47 * Log base2 lookup
48 * l_float32 *makeLogBase2Tab()
49 * l_float32 getLogBase2()
50 *
51 * The image accumulator operations are used when you expect
52 * overflow from 8 bits on intermediate results. For example,
53 * you might want a tophat contrast operator which is
54 * 3*I - opening(I,S) - closing(I,S)
55 * To use these operations, first use the init to generate
56 * a 16 bpp image, use the accumulate to add or subtract 8 bpp
57 * images from that, or the multiply constant to multiply
58 * by a small constant (much less than 256 -- we don't want
59 * overflow from the 16 bit images!), and when you're finished
60 * use final to bring the result back to 8 bpp, clipped
61 * if necessary. There is also a divide function, which
62 * can be used to divide one image by another, scaling the
63 * result for maximum dynamic range, and giving back the
64 * 8 bpp result.
65 *
66 * A simpler interface to the arithmetic operations is
67 * provided in pixacc.c.
68 */
69
70
71 #include <stdio.h>
72 #include <stdlib.h>
73 #include <string.h>
74 #include <math.h>
75 #include "allheaders.h"
76
77
78 /*-------------------------------------------------------------*
79 * One-image grayscale arithmetic operations *
80 *-------------------------------------------------------------*/
81 /*!
82 * pixAddConstantGray()
83 *
84 * Input: pixs (8, 16 or 32 bpp)
85 * val (amount to add to each pixel)
86 * Return: 0 if OK, 1 on error
87 *
88 * Notes:
89 * (1) In-place operation.
90 * (2) No clipping for 32 bpp.
91 * (3) For 8 and 16 bpp, if val > 0 the result is clipped
92 * to 0xff and 0xffff, rsp.
93 * (4) For 8 and 16 bpp, if val < 0 the result is clipped to 0.
94 */
95 l_int32
pixAddConstantGray(PIX * pixs,l_int32 val)96 pixAddConstantGray(PIX *pixs,
97 l_int32 val)
98 {
99 l_int32 w, h, d, wpl;
100 l_uint32 *data;
101
102 PROCNAME("pixAddConstantGray");
103
104 if (!pixs)
105 return ERROR_INT("pixs not defined", procName, 1);
106 pixGetDimensions(pixs, &w, &h, &d);
107 if (d != 8 && d != 16 && d != 32)
108 return ERROR_INT("pixs not 8, 16 or 32 bpp", procName, 1);
109
110 data = pixGetData(pixs);
111 wpl = pixGetWpl(pixs);
112 addConstantGrayLow(data, w, h, d, wpl, val);
113
114 return 0;
115 }
116
117
118 /*!
119 * pixMultConstantGray()
120 *
121 * Input: pixs (8, 16 or 32 bpp)
122 * val (>= 0.0; amount to multiply by each pixel)
123 * Return: 0 if OK, 1 on error
124 *
125 * Notes:
126 * (1) In-place operation; val must be >= 0.
127 * (2) No clipping for 32 bpp.
128 * (3) For 8 and 16 bpp, the result is clipped to 0xff and 0xffff, rsp.
129 */
130 l_int32
pixMultConstantGray(PIX * pixs,l_float32 val)131 pixMultConstantGray(PIX *pixs,
132 l_float32 val)
133 {
134 l_int32 w, h, d, wpl;
135 l_uint32 *data;
136
137 PROCNAME("pixMultConstantGray");
138
139 if (!pixs)
140 return ERROR_INT("pixs not defined", procName, 1);
141 pixGetDimensions(pixs, &w, &h, &d);
142 if (d != 8 && d != 16 && d != 32)
143 return ERROR_INT("pixs not 8, 16 or 32 bpp", procName, 1);
144 if (val < 0.0)
145 return ERROR_INT("val < 0.0", procName, 1);
146
147 data = pixGetData(pixs);
148 wpl = pixGetWpl(pixs);
149 multConstantGrayLow(data, w, h, d, wpl, val);
150
151 return 0;
152 }
153
154
155 /*-------------------------------------------------------------*
156 * Two-image grayscale arithmetic ops *
157 *-------------------------------------------------------------*/
158 /*!
159 * pixAddGray()
160 *
161 * Input: pixd (<optional>; this can be null, equal to pixs1, or
162 * different from pixs1)
163 * pixs1 (can be == to pixd)
164 * pixs2
165 * Return: pixd always
166 *
167 * Notes:
168 * (1) Arithmetic addition of two 8, 16 or 32 bpp images.
169 * (2) For 8 and 16 bpp, we do explicit clipping to 0xff and 0xffff,
170 * respectively.
171 * (3) Alignment is to UL corner.
172 * (4) There are 3 cases. The result can go to a new dest,
173 * in-place to pixs1, or to an existing input dest:
174 * * pixd == null: (src1 + src2) --> new pixd
175 * * pixd == pixs1: (src1 + src2) --> src1 (in-place)
176 * * pixd != pixs1: (src1 + src2) --> input pixd
177 * (5) pixs2 must be different from both pixd and pixs1.
178 */
179 PIX *
pixAddGray(PIX * pixd,PIX * pixs1,PIX * pixs2)180 pixAddGray(PIX *pixd,
181 PIX *pixs1,
182 PIX *pixs2)
183 {
184 l_int32 d, ws, hs, w, h, wpls, wpld;
185 l_uint32 *datas, *datad;
186
187 PROCNAME("pixAddGray");
188
189 if (!pixs1)
190 return (PIX *)ERROR_PTR("pixs1 not defined", procName, pixd);
191 if (!pixs2)
192 return (PIX *)ERROR_PTR("pixs2 not defined", procName, pixd);
193 if (pixs2 == pixs1)
194 return (PIX *)ERROR_PTR("pixs2 and pixs1 must differ", procName, pixd);
195 if (pixs2 == pixd)
196 return (PIX *)ERROR_PTR("pixs2 and pixd must differ", procName, pixd);
197 d = pixGetDepth(pixs1);
198 if (d != 8 && d != 16 && d != 32)
199 return (PIX *)ERROR_PTR("pix are not 8, 16 or 32 bpp", procName, pixd);
200 if (pixGetDepth(pixs2) != d)
201 return (PIX *)ERROR_PTR("depths differ (pixs1, pixs2)", procName, pixd);
202 if (pixd && (pixGetDepth(pixd) != d))
203 return (PIX *)ERROR_PTR("depths differ (pixs1, pixd)", procName, pixd);
204
205 if (!pixSizesEqual(pixs1, pixs2))
206 L_WARNING("pixs1 and pixs2 not equal in size", procName);
207 if (pixd && !pixSizesEqual(pixs1, pixd))
208 L_WARNING("pixs1 and pixd not equal in size", procName);
209
210 if (pixs1 != pixd)
211 pixd = pixCopy(pixd, pixs1);
212
213 /* pixd + pixs2 ==> pixd */
214 datas = pixGetData(pixs2);
215 datad = pixGetData(pixd);
216 wpls = pixGetWpl(pixs2);
217 wpld = pixGetWpl(pixd);
218 pixGetDimensions(pixs2, &ws, &hs, NULL);
219 pixGetDimensions(pixd, &w, &h, NULL);
220 w = L_MIN(ws, w);
221 h = L_MIN(hs, h);
222 addGrayLow(datad, w, h, d, wpld, datas, wpls);
223
224 return pixd;
225 }
226
227
228 /*!
229 * pixSubtractGray()
230 *
231 * Input: pixd (<optional>; this can be null, equal to pixs1, or
232 * different from pixs1)
233 * pixs1 (can be == to pixd)
234 * pixs2
235 * Return: pixd always
236 *
237 * Notes:
238 * (1) Arithmetic subtraction of two 8, 16 or 32 bpp images.
239 * (2) Source pixs2 is always subtracted from source pixs1.
240 * (3) Do explicit clipping to 0.
241 * (4) Alignment is to UL corner.
242 * (5) There are 3 cases. The result can go to a new dest,
243 * in-place to pixs1, or to an existing input dest:
244 * (a) pixd == null (src1 - src2) --> new pixd
245 * (b) pixd == pixs1 (src1 - src2) --> src1 (in-place)
246 * (d) pixd != pixs1 (src1 - src2) --> input pixd
247 * (6) pixs2 must be different from both pixd and pixs1.
248 */
249 PIX *
pixSubtractGray(PIX * pixd,PIX * pixs1,PIX * pixs2)250 pixSubtractGray(PIX *pixd,
251 PIX *pixs1,
252 PIX *pixs2)
253 {
254 l_int32 w, h, ws, hs, d, wpls, wpld;
255 l_uint32 *datas, *datad;
256
257 PROCNAME("pixSubtractGray");
258
259 if (!pixs1)
260 return (PIX *)ERROR_PTR("pixs1 not defined", procName, pixd);
261 if (!pixs2)
262 return (PIX *)ERROR_PTR("pixs2 not defined", procName, pixd);
263 if (pixs2 == pixs1)
264 return (PIX *)ERROR_PTR("pixs2 and pixs1 must differ", procName, pixd);
265 if (pixs2 == pixd)
266 return (PIX *)ERROR_PTR("pixs2 and pixd must differ", procName, pixd);
267 d = pixGetDepth(pixs1);
268 if (d != 8 && d != 16 && d != 32)
269 return (PIX *)ERROR_PTR("pix are not 8, 16 or 32 bpp", procName, pixd);
270 if (pixGetDepth(pixs2) != d)
271 return (PIX *)ERROR_PTR("depths differ (pixs1, pixs2)", procName, pixd);
272 if (pixd && (pixGetDepth(pixd) != d))
273 return (PIX *)ERROR_PTR("depths differ (pixs1, pixd)", procName, pixd);
274
275 if (!pixSizesEqual(pixs1, pixs2))
276 L_WARNING("pixs1 and pixs2 not equal in size", procName);
277 if (pixd && !pixSizesEqual(pixs1, pixd))
278 L_WARNING("pixs1 and pixd not equal in size", procName);
279
280 if (pixs1 != pixd)
281 pixd = pixCopy(pixd, pixs1);
282
283 /* pixd - pixs2 ==> pixd */
284 datas = pixGetData(pixs2);
285 datad = pixGetData(pixd);
286 wpls = pixGetWpl(pixs2);
287 wpld = pixGetWpl(pixd);
288 pixGetDimensions(pixs2, &ws, &hs, NULL);
289 pixGetDimensions(pixd, &w, &h, NULL);
290 w = L_MIN(ws, w);
291 h = L_MIN(hs, h);
292 subtractGrayLow(datad, w, h, d, wpld, datas, wpls);
293
294 return pixd;
295 }
296
297
298 /*-------------------------------------------------------------*
299 * Grayscale threshold operation *
300 *-------------------------------------------------------------*/
301 /*!
302 * pixThresholdToValue()
303 *
304 * Input: pixd (<optional>; if not null, must be equal to pixs)
305 * pixs (8, 16, 32 bpp)
306 * threshval
307 * setval
308 * Return: pixd always
309 *
310 * Notes:
311 * - operation can be in-place (pixs == pixd) or to a new pixd
312 * - if setval > threshval, sets pixels with a value >= threshval to setval
313 * - if setval < threshval, sets pixels with a value <= threshval to setval
314 * - if setval == threshval, no-op
315 */
316 PIX *
pixThresholdToValue(PIX * pixd,PIX * pixs,l_int32 threshval,l_int32 setval)317 pixThresholdToValue(PIX *pixd,
318 PIX *pixs,
319 l_int32 threshval,
320 l_int32 setval)
321 {
322 l_int32 w, h, d, wpld;
323 l_uint32 *datad;
324
325 PROCNAME("pixThresholdToValue");
326
327 if (!pixs)
328 return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
329 d = pixGetDepth(pixs);
330 if (d != 8 && d != 16 && d != 32)
331 return (PIX *)ERROR_PTR("pixs not 8, 16 or 32 bpp", procName, pixd);
332 if (pixd && (pixs != pixd))
333 return (PIX *)ERROR_PTR("pixd exists and is not pixs", procName, pixd);
334 if (threshval < 0 || setval < 0)
335 return (PIX *)ERROR_PTR("threshval & setval not < 0", procName, pixd);
336 if (d == 8 && setval > 255)
337 return (PIX *)ERROR_PTR("setval > 255 for 8 bpp", procName, pixd);
338 if (d == 16 && setval > 0xffff)
339 return (PIX *)ERROR_PTR("setval > 0xffff for 16 bpp", procName, pixd);
340
341 if (!pixd)
342 pixd = pixCopy(NULL, pixs);
343 if (setval == threshval) {
344 L_WARNING("setval == threshval; no operation", procName);
345 return pixd;
346 }
347
348 datad = pixGetData(pixd);
349 pixGetDimensions(pixd, &w, &h, NULL);
350 wpld = pixGetWpl(pixd);
351
352 thresholdToValueLow(datad, w, h, d, wpld, threshval, setval);
353 return pixd;
354 }
355
356
357
358 /*-------------------------------------------------------------*
359 * Image accumulator arithmetic operations *
360 *-------------------------------------------------------------*/
361 /*!
362 * pixInitAccumulate()
363 *
364 * Input: w, h (of accumulate array)
365 * offset (initialize the 32 bpp to have this
366 * value; not more than 0x40000000)
367 * Return: pixd (32 bpp), or null on error
368 *
369 * Notes:
370 * (1) The offset must be >= 0.
371 * (2) The offset is used so that we can do arithmetic
372 * with negative number results on l_uint32 data; it
373 * prevents the l_uint32 data from going negative.
374 * (3) Because we use l_int32 intermediate data results,
375 * these should never exceed the max of l_int32 (0x7fffffff).
376 * We do not permit the offset to be above 0x40000000,
377 * which is half way between 0 and the max of l_int32.
378 * (4) The same offset should be used for initialization,
379 * multiplication by a constant, and final extraction!
380 * (5) If you're only adding positive values, offset can be 0.
381 */
382 PIX *
pixInitAccumulate(l_int32 w,l_int32 h,l_uint32 offset)383 pixInitAccumulate(l_int32 w,
384 l_int32 h,
385 l_uint32 offset)
386 {
387 PIX *pixd;
388
389 PROCNAME("pixInitAccumulate");
390
391 if ((pixd = pixCreate(w, h, 32)) == NULL)
392 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
393 if (offset > 0x40000000)
394 offset = 0x40000000;
395 pixSetAllArbitrary(pixd, offset);
396 return pixd;
397 }
398
399
400 /*!
401 * pixFinalAccumulate()
402 *
403 * Input: pixs (32 bpp)
404 * offset (same as used for initialization)
405 * depth (8, 16 or 32 bpp, of destination)
406 * Return: pixd (8, 16 or 32 bpp), or null on error
407 *
408 * Notes:
409 * (1) The offset must be >= 0 and should not exceed 0x40000000.
410 * (2) The offset is subtracted from the src 32 bpp image
411 * (3) For 8 bpp dest, the result is clipped to [0, 0xff]
412 * (4) For 16 bpp dest, the result is clipped to [0, 0xffff]
413 */
414 PIX *
pixFinalAccumulate(PIX * pixs,l_uint32 offset,l_int32 depth)415 pixFinalAccumulate(PIX *pixs,
416 l_uint32 offset,
417 l_int32 depth)
418 {
419 l_int32 w, h, wpls, wpld;
420 l_uint32 *datas, *datad;
421 PIX *pixd;
422
423 PROCNAME("pixFinalAccumulate");
424
425 if (!pixs)
426 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
427 if (pixGetDepth(pixs) != 32)
428 return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
429 if (depth != 8 && depth != 16 && depth != 32)
430 return (PIX *)ERROR_PTR("dest depth not 8, 16, 32 bpp", procName, NULL);
431 if (offset > 0x40000000)
432 offset = 0x40000000;
433
434 pixGetDimensions(pixs, &w, &h, NULL);
435 if ((pixd = pixCreate(w, h, depth)) == NULL)
436 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
437 pixCopyResolution(pixd, pixs); /* but how did pixs get it initially? */
438 datas = pixGetData(pixs);
439 datad = pixGetData(pixd);
440 wpls = pixGetWpl(pixs);
441 wpld = pixGetWpl(pixd);
442
443 finalAccumulateLow(datad, w, h, depth, wpld, datas, wpls, offset);
444 return pixd;
445 }
446
447
448 /*!
449 * pixFinalAccumulateThreshold()
450 *
451 * Input: pixs (32 bpp)
452 * offset (same as used for initialization)
453 * threshold (values less than this are set in the destination)
454 * Return: pixd (1 bpp), or null on error
455 *
456 * Notes:
457 * (1) The offset must be >= 0 and should not exceed 0x40000000.
458 * (2) The offset is subtracted from the src 32 bpp image
459 */
460 PIX *
pixFinalAccumulateThreshold(PIX * pixs,l_uint32 offset,l_uint32 threshold)461 pixFinalAccumulateThreshold(PIX *pixs,
462 l_uint32 offset,
463 l_uint32 threshold)
464 {
465 l_int32 w, h, wpls, wpld;
466 l_uint32 *datas, *datad;
467 PIX *pixd;
468
469 PROCNAME("pixFinalAccumulateThreshold");
470
471 if (!pixs)
472 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
473 if (pixGetDepth(pixs) != 32)
474 return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
475 if (offset > 0x40000000)
476 offset = 0x40000000;
477
478 pixGetDimensions(pixs, &w, &h, NULL);
479 if ((pixd = pixCreate(w, h, 1)) == NULL)
480 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
481 pixCopyResolution(pixd, pixs); /* but how did pixs get it initially? */
482 datas = pixGetData(pixs);
483 datad = pixGetData(pixd);
484 wpls = pixGetWpl(pixs);
485 wpld = pixGetWpl(pixd);
486
487 finalAccumulateThreshLow(datad, w, h, wpld, datas, wpls, offset, threshold);
488 return pixd;
489 }
490
491
492 /*!
493 * pixAccumulate()
494 *
495 * Input: pixd (32 bpp)
496 * pixs (1, 8, 16 or 32 bpp)
497 * op (L_ARITH_ADD or L_ARITH_SUBTRACT)
498 * Return: 0 if OK; 1 on error
499 *
500 * Notes:
501 * (1) This adds or subtracts each pixs value from pixd.
502 * (2) This clips to the minimum of pixs and pixd, so they
503 * do not need to be the same size.
504 * (3) The alignment is to the origin (UL corner) of pixs & pixd.
505 */
506 l_int32
pixAccumulate(PIX * pixd,PIX * pixs,l_int32 op)507 pixAccumulate(PIX *pixd,
508 PIX *pixs,
509 l_int32 op)
510 {
511 l_int32 w, h, d, wd, hd, wpls, wpld;
512 l_uint32 *datas, *datad;
513
514 PROCNAME("pixAccumulate");
515
516 if (!pixd || (pixGetDepth(pixd) != 32))
517 return ERROR_INT("pixd not defined or not 32 bpp", procName, 1);
518 if (!pixs)
519 return ERROR_INT("pixs not defined", procName, 1);
520 d = pixGetDepth(pixs);
521 if (d != 1 && d != 8 && d != 16 && d != 32)
522 return ERROR_INT("pixs not 1, 8, 16 or 32 bpp", procName, 1);
523 if (op != L_ARITH_ADD && op != L_ARITH_SUBTRACT)
524 return ERROR_INT("op must be in {L_ARITH_ADD, L_ARITH_SUBTRACT}",
525 procName, 1);
526
527 datas = pixGetData(pixs);
528 datad = pixGetData(pixd);
529 wpls = pixGetWpl(pixs);
530 wpld = pixGetWpl(pixd);
531 pixGetDimensions(pixs, &w, &h, NULL);
532 pixGetDimensions(pixd, &wd, &hd, NULL);
533 w = L_MIN(w, wd);
534 h = L_MIN(h, hd);
535
536 accumulateLow(datad, w, h, wpld, datas, d, wpls, op);
537 return 0;
538 }
539
540
541 /*!
542 * pixMultConstAccumulate()
543 *
544 * Input: pixs (32 bpp)
545 * factor
546 * offset (same as used for initialization)
547 * Return: 0 if OK; 1 on error
548 *
549 * Notes:
550 * (1) The offset must be >= 0 and should not exceed 0x40000000.
551 * (2) This multiplies each pixel, relative to offset, by the input factor
552 * (3) The result is returned with the offset back in place.
553 */
554 l_int32
pixMultConstAccumulate(PIX * pixs,l_float32 factor,l_uint32 offset)555 pixMultConstAccumulate(PIX *pixs,
556 l_float32 factor,
557 l_uint32 offset)
558 {
559 l_int32 w, h, wpl;
560 l_uint32 *data;
561
562 PROCNAME("pixMultConstAccumulate");
563
564 if (!pixs)
565 return ERROR_INT("pixs not defined", procName, 1);
566 if (pixGetDepth(pixs) != 32)
567 return ERROR_INT("pixs not 32 bpp", procName, 1);
568 if (offset > 0x40000000)
569 offset = 0x40000000;
570
571 pixGetDimensions(pixs, &w, &h, NULL);
572 data = pixGetData(pixs);
573 wpl = pixGetWpl(pixs);
574
575 multConstAccumulateLow(data, w, h, wpl, factor, offset);
576 return 0;
577 }
578
579
580 /*-----------------------------------------------------------------------*
581 * Absolute value of difference *
582 *-----------------------------------------------------------------------*/
583 /*!
584 * pixAbsDifference()
585 *
586 * Input: pixs1, pixs2 (both either 8 or 16 bpp gray, or 32 bpp RGB)
587 * Return: pixd, or null on error
588 *
589 * Notes:
590 * (1) The depth of pixs1 and pixs2 must be equal.
591 * (2) Clips computation to the min size, aligning the UL corners
592 * (3) For 8 and 16 bpp, assumes one gray component.
593 * (4) For 32 bpp, assumes 3 color components, and ignores the
594 * LSB of each word (the alpha channel)
595 * (5) Computes the absolute value of the difference between
596 * each component value.
597 */
598 PIX *
pixAbsDifference(PIX * pixs1,PIX * pixs2)599 pixAbsDifference(PIX *pixs1,
600 PIX *pixs2)
601 {
602 l_int32 w, h, w2, h2, d, wpls, wpld;
603 l_uint32 *datas1, *datas2, *datad;
604 PIX *pixd;
605
606 PROCNAME("pixAbsDifference");
607
608 if (!pixs1)
609 return (PIX *)ERROR_PTR("pixs1 not defined", procName, NULL);
610 if (!pixs2)
611 return (PIX *)ERROR_PTR("pixs2 not defined", procName, NULL);
612 d = pixGetDepth(pixs1);
613 if (d != pixGetDepth(pixs2))
614 return (PIX *)ERROR_PTR("src1 and src2 depths unequal", procName, NULL);
615 if (d != 8 && d != 16 && d != 32)
616 return (PIX *)ERROR_PTR("depths not in {8, 16, 32}", procName, NULL);
617
618 pixGetDimensions(pixs1, &w, &h, NULL);
619 pixGetDimensions(pixs2, &w2, &h2, NULL);
620 w = L_MIN(w, w2);
621 h = L_MIN(h, h2);
622 if ((pixd = pixCreate(w, h, d)) == NULL)
623 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
624 pixCopyResolution(pixd, pixs1);
625 datas1 = pixGetData(pixs1);
626 datas2 = pixGetData(pixs2);
627 datad = pixGetData(pixd);
628 wpls = pixGetWpl(pixs1);
629 wpld = pixGetWpl(pixd);
630
631 absDifferenceLow(datad, w, h, wpld, datas1, datas2, d, wpls);
632
633 return pixd;
634 }
635
636
637 /*-----------------------------------------------------------------------*
638 * Two-image min and max operations (8 and 16 bpp) *
639 *-----------------------------------------------------------------------*/
640 /*!
641 * pixMinOrMax()
642 *
643 * Input: pixd (<optional> destination: this can be null,
644 * equal to pixs1, or different from pixs1)
645 * pixs1 (can be == to pixd)
646 * pixs2
647 * type (L_CHOOSE_MIN, L_CHOOSE_MAX)
648 * Return: pixd always
649 *
650 * Notes:
651 * (1) This gives the min or max of two images.
652 * (2) The depth can be 8 or 16 bpp.
653 * (3) There are 3 cases:
654 * - if pixd == null, Min(src1, src2) --> new pixd
655 * - if pixd == pixs1, Min(src1, src2) --> src1 (in-place)
656 * - if pixd != pixs1, Min(src1, src2) --> input pixd
657 */
658 PIX *
pixMinOrMax(PIX * pixd,PIX * pixs1,PIX * pixs2,l_int32 type)659 pixMinOrMax(PIX *pixd,
660 PIX *pixs1,
661 PIX *pixs2,
662 l_int32 type)
663 {
664 l_int32 d, ws, hs, w, h, wpls, wpld, i, j;
665 l_int32 vals, vald, val;
666 l_uint32 *datas, *datad, *lines, *lined;
667
668 PROCNAME("pixMinOrMax");
669
670 if (!pixs1)
671 return (PIX *)ERROR_PTR("pixs1 not defined", procName, pixd);
672 if (!pixs2)
673 return (PIX *)ERROR_PTR("pixs2 not defined", procName, pixd);
674 if (pixs1 == pixs2)
675 return (PIX *)ERROR_PTR("pixs1 and pixs2 must differ", procName, pixd);
676 if (type != L_CHOOSE_MIN && type != L_CHOOSE_MAX)
677 return (PIX *)ERROR_PTR("invalid type", procName, pixd);
678 d = pixGetDepth(pixs1);
679 if (pixGetDepth(pixs2) != d)
680 return (PIX *)ERROR_PTR("depths unequal", procName, pixd);
681 if (d != 8 && d != 16)
682 return (PIX *)ERROR_PTR("depth not 8 or 16 bpp", procName, pixd);
683
684 if (pixs1 != pixd)
685 pixd = pixCopy(pixd, pixs1);
686
687 pixGetDimensions(pixs2, &ws, &hs, NULL);
688 pixGetDimensions(pixd, &w, &h, NULL);
689 w = L_MIN(w, ws);
690 h = L_MIN(h, hs);
691 datas = pixGetData(pixs2);
692 datad = pixGetData(pixd);
693 wpls = pixGetWpl(pixs2);
694 wpld = pixGetWpl(pixd);
695 for (i = 0; i < h; i++) {
696 lines = datas + i * wpls;
697 lined = datad + i * wpld;
698 if (d == 8) {
699 if (type == L_CHOOSE_MIN) {
700 for (j = 0; j < w; j++) {
701 vals = GET_DATA_BYTE(lines, j);
702 vald = GET_DATA_BYTE(lined, j);
703 val = L_MIN(vals, vald);
704 SET_DATA_BYTE(lined, j, val);
705 }
706 } else { /* type == L_CHOOSE_MAX */
707 for (j = 0; j < w; j++) {
708 vals = GET_DATA_BYTE(lines, j);
709 vald = GET_DATA_BYTE(lined, j);
710 val = L_MAX(vals, vald);
711 SET_DATA_BYTE(lined, j, val);
712 }
713 }
714 } else { /* d == 16 */
715 if (type == L_CHOOSE_MIN) {
716 for (j = 0; j < w; j++) {
717 vals = GET_DATA_TWO_BYTES(lines, j);
718 vald = GET_DATA_TWO_BYTES(lined, j);
719 val = L_MIN(vals, vald);
720 SET_DATA_TWO_BYTES(lined, j, val);
721 }
722 } else { /* type == L_CHOOSE_MAX */
723 for (j = 0; j < w; j++) {
724 vals = GET_DATA_TWO_BYTES(lines, j);
725 vald = GET_DATA_TWO_BYTES(lined, j);
726 val = L_MAX(vals, vald);
727 SET_DATA_TWO_BYTES(lined, j, val);
728 }
729 }
730 }
731 }
732
733 return pixd;
734 }
735
736
737 /*-----------------------------------------------------------------------*
738 * Scale for maximum dynamic range in 8 bpp image *
739 *-----------------------------------------------------------------------*/
740 /*!
741 * pixMaxDynamicRange()
742 *
743 * Input: pixs (4, 8, 16 or 32 bpp source)
744 * type (L_LINEAR_SCALE or L_LOG_SCALE)
745 * Return: pixd (8 bpp), or null on error
746 *
747 * Notes:
748 * (1) Scales pixel values to fit maximally within the dest 8 bpp pixd
749 * (2) Uses a LUT for log scaling
750 */
751 PIX *
pixMaxDynamicRange(PIX * pixs,l_int32 type)752 pixMaxDynamicRange(PIX *pixs,
753 l_int32 type)
754 {
755 l_uint8 dval;
756 l_int32 i, j, w, h, d, wpls, wpld, max, sval;
757 l_uint32 *datas, *datad;
758 l_uint32 word;
759 l_uint32 *lines, *lined;
760 l_float32 factor;
761 l_float32 *tab;
762 PIX *pixd;
763
764 PROCNAME("pixMaxDynamicRange");
765
766 if (!pixs)
767 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
768 d = pixGetDepth(pixs);
769 if (d != 4 && d != 8 && d != 16 && d != 32)
770 return (PIX *)ERROR_PTR("pixs not in {4,8,16,32} bpp", procName, NULL);
771 if (type != L_LINEAR_SCALE && type != L_LOG_SCALE)
772 return (PIX *)ERROR_PTR("invalid type", procName, NULL);
773
774 pixGetDimensions(pixs, &w, &h, NULL);
775 if ((pixd = pixCreate(w, h, 8)) == NULL)
776 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
777 pixCopyResolution(pixd, pixs);
778 datas = pixGetData(pixs);
779 datad = pixGetData(pixd);
780 wpls = pixGetWpl(pixs);
781 wpld = pixGetWpl(pixd);
782
783 /* Get max */
784 max = 0;
785 for (i = 0; i < h; i++) {
786 lines = datas + i * wpls;
787 for (j = 0; j < wpls; j++) {
788 word = *(lines + j);
789 if (d == 4) {
790 max = L_MAX(max, word >> 28);
791 max = L_MAX(max, (word >> 24) & 0xf);
792 max = L_MAX(max, (word >> 20) & 0xf);
793 max = L_MAX(max, (word >> 16) & 0xf);
794 max = L_MAX(max, (word >> 12) & 0xf);
795 max = L_MAX(max, (word >> 8) & 0xf);
796 max = L_MAX(max, (word >> 4) & 0xf);
797 max = L_MAX(max, word & 0xf);
798 } else if (d == 8) {
799 max = L_MAX(max, word >> 24);
800 max = L_MAX(max, (word >> 16) & 0xff);
801 max = L_MAX(max, (word >> 8) & 0xff);
802 max = L_MAX(max, word & 0xff);
803 } else if (d == 16) {
804 max = L_MAX(max, word >> 16);
805 max = L_MAX(max, word & 0xffff);
806 } else { /* d == 32 */
807 max = L_MAX(max, word);
808 }
809 }
810 }
811
812 /* Map to the full dynamic range of 8 bpp output */
813 if (d == 4) {
814 if (type == L_LINEAR_SCALE) {
815 factor = 255. / (l_float32)max;
816 for (i = 0; i < h; i++) {
817 lines = datas + i * wpls;
818 lined = datad + i * wpld;
819 for (j = 0; j < w; j++) {
820 sval = GET_DATA_QBIT(lines, j);
821 dval = (l_uint8)(factor * (l_float32)sval + 0.5);
822 SET_DATA_QBIT(lined, j, dval);
823 }
824 }
825 } else { /* type == L_LOG_SCALE) */
826 tab = makeLogBase2Tab();
827 factor = 255. / getLogBase2(max, tab);
828 for (i = 0; i < h; i++) {
829 lines = datas + i * wpls;
830 lined = datad + i * wpld;
831 for (j = 0; j < w; j++) {
832 sval = GET_DATA_QBIT(lines, j);
833 dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
834 SET_DATA_BYTE(lined, j, dval);
835 }
836 }
837 FREE(tab);
838 }
839 } else if (d == 8) {
840 if (type == L_LINEAR_SCALE) {
841 factor = 255. / (l_float32)max;
842 for (i = 0; i < h; i++) {
843 lines = datas + i * wpls;
844 lined = datad + i * wpld;
845 for (j = 0; j < w; j++) {
846 sval = GET_DATA_BYTE(lines, j);
847 dval = (l_uint8)(factor * (l_float32)sval + 0.5);
848 SET_DATA_BYTE(lined, j, dval);
849 }
850 }
851 } else { /* type == L_LOG_SCALE) */
852 tab = makeLogBase2Tab();
853 factor = 255. / getLogBase2(max, tab);
854 for (i = 0; i < h; i++) {
855 lines = datas + i * wpls;
856 lined = datad + i * wpld;
857 for (j = 0; j < w; j++) {
858 sval = GET_DATA_BYTE(lines, j);
859 dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
860 SET_DATA_BYTE(lined, j, dval);
861 }
862 }
863 FREE(tab);
864 }
865 } else if (d == 16) {
866 if (type == L_LINEAR_SCALE) {
867 factor = 255. / (l_float32)max;
868 for (i = 0; i < h; i++) {
869 lines = datas + i * wpls;
870 lined = datad + i * wpld;
871 for (j = 0; j < w; j++) {
872 sval = GET_DATA_TWO_BYTES(lines, j);
873 dval = (l_uint8)(factor * (l_float32)sval + 0.5);
874 SET_DATA_BYTE(lined, j, dval);
875 }
876 }
877 } else { /* type == L_LOG_SCALE) */
878 tab = makeLogBase2Tab();
879 factor = 255. / getLogBase2(max, tab);
880 for (i = 0; i < h; i++) {
881 lines = datas + i * wpls;
882 lined = datad + i * wpld;
883 for (j = 0; j < w; j++) {
884 sval = GET_DATA_TWO_BYTES(lines, j);
885 dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
886 SET_DATA_BYTE(lined, j, dval);
887 }
888 }
889 FREE(tab);
890 }
891 } else { /* d == 32 */
892 if (type == L_LINEAR_SCALE) {
893 factor = 255. / (l_float32)max;
894 for (i = 0; i < h; i++) {
895 lines = datas + i * wpls;
896 lined = datad + i * wpld;
897 for (j = 0; j < w; j++) {
898 sval = lines[j];
899 dval = (l_uint8)(factor * (l_float32)sval + 0.5);
900 SET_DATA_BYTE(lined, j, dval);
901 }
902 }
903 } else { /* type == L_LOG_SCALE) */
904 tab = makeLogBase2Tab();
905 factor = 255. / getLogBase2(max, tab);
906 for (i = 0; i < h; i++) {
907 lines = datas + i * wpls;
908 lined = datad + i * wpld;
909 for (j = 0; j < w; j++) {
910 sval = lines[j];
911 dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
912 SET_DATA_BYTE(lined, j, dval);
913 }
914 }
915 FREE(tab);
916 }
917 }
918
919 return pixd;
920 }
921
922
923 /*-----------------------------------------------------------------------*
924 * Log base2 lookup *
925 *-----------------------------------------------------------------------*/
926 /*
927 * makeLogBase2Tab()
928 *
929 * Input: void
930 * Return: table (giving the log[base 2] of val)
931 */
932 l_float32 *
makeLogBase2Tab(void)933 makeLogBase2Tab(void)
934 {
935 l_int32 i;
936 l_float32 log2;
937 l_float32 *tab;
938
939 PROCNAME("makeLogBase2Tab");
940
941 if ((tab = (l_float32 *)CALLOC(256, sizeof(l_float32))) == NULL)
942 return (l_float32 *)ERROR_PTR("tab not made", procName, NULL);
943
944 log2 = (l_float32)log((l_float32)2);
945 for (i = 0; i < 256; i++)
946 tab[i] = (l_float32)log((l_float32)i) / log2;
947
948 return tab;
949 }
950
951
952 /*
953 * getLogBase2()
954 *
955 * Input: val
956 * logtab (256-entry table of logs)
957 * Return: logdist, or 0 on error
958 */
959 l_float32
getLogBase2(l_int32 val,l_float32 * logtab)960 getLogBase2(l_int32 val,
961 l_float32 *logtab)
962 {
963 PROCNAME("getLogBase2");
964
965 if (!logtab)
966 return ERROR_INT("logtab not defined", procName, 0);
967
968 if (val < 0x100)
969 return logtab[val];
970 else if (val < 0x10000)
971 return 8.0 + logtab[val >> 8];
972 else if (val < 0x1000000)
973 return 16.0 + logtab[val >> 16];
974 else
975 return 24.0 + logtab[val >> 24];
976 }
977
978