Mercurial > hgrepos > Python2 > PyMuPDF

comparison mupdf-source/thirdparty/leptonica/src/kernel.c @ 2:b50eed0cc0ef upstream

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ADD: MuPDF v1.26.7: the MuPDF source as downloaded by a default build of PyMuPDF 1.26.4. The directory name has changed: no version number in the expanded directory now.
author Franz Glasner <fzglas.hg@dom66.de>
date Mon, 15 Sep 2025 11:43:07 +0200
parents
children
comparison
equal deleted inserted replaced
1:1d09e1dec1d9 2:b50eed0cc0ef
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 /*!
29 * \file kernel.c
30 * <pre>
31 *
32 * Basic operations on kernels for image convolution
33 *
34 * Create/destroy/copy
35 * L_KERNEL *kernelCreate()
36 * void kernelDestroy()
37 * L_KERNEL *kernelCopy()
38 *
39 * Accessors:
40 * l_int32 kernelGetElement()
41 * l_int32 kernelSetElement()
42 * l_int32 kernelGetParameters()
43 * l_int32 kernelSetOrigin()
44 * l_int32 kernelGetSum()
45 * l_int32 kernelGetMinMax()
46 *
47 * Normalize/invert
48 * L_KERNEL *kernelNormalize()
49 * L_KERNEL *kernelInvert()
50 *
51 * Helper function
52 * l_float32 **create2dFloatArray()
53 *
54 * Serialized I/O
55 * L_KERNEL *kernelRead()
56 * L_KERNEL *kernelReadStream()
57 * l_int32 kernelWrite()
58 * l_int32 kernelWriteStream()
59 *
60 * Making a kernel from a compiled string
61 * L_KERNEL *kernelCreateFromString()
62 *
63 * Making a kernel from a simple file format
64 * L_KERNEL *kernelCreateFromFile()
65 *
66 * Making a kernel from a Pix
67 * L_KERNEL *kernelCreateFromPix()
68 *
69 * Display a kernel in a pix
70 * PIX *kernelDisplayInPix()
71 *
72 * Parse string to extract numbers
73 * NUMA *parseStringForNumbers()
74 *
75 * Simple parametric kernels
76 * L_KERNEL *makeFlatKernel()
77 * L_KERNEL *makeGaussianKernel()
78 * L_KERNEL *makeGaussianKernelSep()
79 * L_KERNEL *makeDoGKernel()
80 * </pre>
81 */
82
83 #ifdef HAVE_CONFIG_H
84 #include <config_auto.h>
85 #endif /* HAVE_CONFIG_H */
86
87 #include <string.h>
88 #include <math.h>
89 #include "allheaders.h"
90
91 /* Array size must be > 0 and not larger than this */
92 static const l_uint32 MaxArraySize = 100000;
93
94 /*------------------------------------------------------------------------*
95 * Create / Destroy *
96 *------------------------------------------------------------------------*/
97 /*!
98 * \brief kernelCreate()
99 *
100 * \param[in] height, width
101 * \return kernel, or NULL on error
102 *
103 * <pre>
104 * Notes:
105 * (1) kernelCreate() initializes all values to 0.
106 * (2) After this call, (cy,cx) and nonzero data values must be
107 * assigned.
108 * (2) The number of kernel elements must be less than 2^29.
109 * </pre>
110 */
111 L_KERNEL *
112 kernelCreate(l_int32 height,
113 l_int32 width)
114 {
115 l_uint64 size64;
116 L_KERNEL *kel;
117
118 if (width <= 0)
119 return (L_KERNEL *)ERROR_PTR("width must be > 0", __func__, NULL);
120 if (height <= 0)
121 return (L_KERNEL *)ERROR_PTR("height must be > 0", __func__, NULL);
122
123 /* Avoid overflow in malloc arg */
124 size64 = (l_uint64)width * (l_uint64)height;
125 if (size64 >= (1LL << 29)) {
126 L_ERROR("requested width = %d, height = %d\n", __func__, width, height);
127 return (L_KERNEL *)ERROR_PTR("size >= 2^29", __func__, NULL);
128 }
129
130 kel = (L_KERNEL *)LEPT_CALLOC(1, sizeof(L_KERNEL));
131 kel->sy = height;
132 kel->sx = width;
133 if ((kel->data = create2dFloatArray(height, width)) == NULL) {
134 LEPT_FREE(kel);
135 return (L_KERNEL *)ERROR_PTR("data not allocated", __func__, NULL);
136 }
137 return kel;
138 }
139
140
141 /*!
142 * \brief kernelDestroy()
143 *
144 * \param[in,out] pkel will be set to null before returning
145 * \return void
146 */
147 void
148 kernelDestroy(L_KERNEL **pkel)
149 {
150 l_int32 i;
151 L_KERNEL *kel;
152
153 if (pkel == NULL) {
154 L_WARNING("ptr address is NULL!\n", __func__);
155 return;
156 }
157 if ((kel = *pkel) == NULL)
158 return;
159
160 for (i = 0; i < kel->sy; i++)
161 LEPT_FREE(kel->data[i]);
162 LEPT_FREE(kel->data);
163 LEPT_FREE(kel);
164 *pkel = NULL;
165 }
166
167
168 /*!
169 * \brief kernelCopy()
170 *
171 * \param[in] kels source kernel
172 * \return keld copy of kels, or NULL on error
173 */
174 L_KERNEL *
175 kernelCopy(L_KERNEL *kels)
176 {
177 l_int32 i, j, sx, sy, cx, cy;
178 L_KERNEL *keld;
179
180 if (!kels)
181 return (L_KERNEL *)ERROR_PTR("kels not defined", __func__, NULL);
182
183 kernelGetParameters(kels, &sy, &sx, &cy, &cx);
184 if ((keld = kernelCreate(sy, sx)) == NULL)
185 return (L_KERNEL *)ERROR_PTR("keld not made", __func__, NULL);
186 keld->cy = cy;
187 keld->cx = cx;
188 for (i = 0; i < sy; i++)
189 for (j = 0; j < sx; j++)
190 keld->data[i][j] = kels->data[i][j];
191
192 return keld;
193 }
194
195
196 /*----------------------------------------------------------------------*
197 * Accessors *
198 *----------------------------------------------------------------------*/
199 /*!
200 * \brief kernelGetElement()
201 *
202 * \param[in] kel
203 * \param[in] row
204 * \param[in] col
205 * \param[out] pval
206 * \return 0 if OK; 1 on error
207 */
208 l_ok
209 kernelGetElement(L_KERNEL *kel,
210 l_int32 row,
211 l_int32 col,
212 l_float32 *pval)
213 {
214 if (!pval)
215 return ERROR_INT("&val not defined", __func__, 1);
216 *pval = 0;
217 if (!kel)
218 return ERROR_INT("kernel not defined", __func__, 1);
219 if (row < 0 || row >= kel->sy)
220 return ERROR_INT("kernel row out of bounds", __func__, 1);
221 if (col < 0 || col >= kel->sx)
222 return ERROR_INT("kernel col out of bounds", __func__, 1);
223
224 *pval = kel->data[row][col];
225 return 0;
226 }
227
228
229 /*!
230 * \brief kernelSetElement()
231 *
232 * \param[in] kel kernel
233 * \param[in] row
234 * \param[in] col
235 * \param[in] val
236 * \return 0 if OK; 1 on error
237 */
238 l_ok
239 kernelSetElement(L_KERNEL *kel,
240 l_int32 row,
241 l_int32 col,
242 l_float32 val)
243 {
244 if (!kel)
245 return ERROR_INT("kel not defined", __func__, 1);
246 if (row < 0 || row >= kel->sy)
247 return ERROR_INT("kernel row out of bounds", __func__, 1);
248 if (col < 0 || col >= kel->sx)
249 return ERROR_INT("kernel col out of bounds", __func__, 1);
250
251 kel->data[row][col] = val;
252 return 0;
253 }
254
255
256 /*!
257 * \brief kernelGetParameters()
258 *
259 * \param[in] kel kernel
260 * \param[out] psy, psx, pcy, pcx [optional] each can be null
261 * \return 0 if OK, 1 on error
262 */
263 l_ok
264 kernelGetParameters(L_KERNEL *kel,
265 l_int32 *psy,
266 l_int32 *psx,
267 l_int32 *pcy,
268 l_int32 *pcx)
269 {
270 if (psy) *psy = 0;
271 if (psx) *psx = 0;
272 if (pcy) *pcy = 0;
273 if (pcx) *pcx = 0;
274 if (!kel)
275 return ERROR_INT("kernel not defined", __func__, 1);
276 if (psy) *psy = kel->sy;
277 if (psx) *psx = kel->sx;
278 if (pcy) *pcy = kel->cy;
279 if (pcx) *pcx = kel->cx;
280 return 0;
281 }
282
283
284 /*!
285 * \brief kernelSetOrigin()
286 *
287 * \param[in] kel kernel
288 * \param[in] cy, cx
289 * \return 0 if OK; 1 on error
290 */
291 l_ok
292 kernelSetOrigin(L_KERNEL *kel,
293 l_int32 cy,
294 l_int32 cx)
295 {
296 if (!kel)
297 return ERROR_INT("kel not defined", __func__, 1);
298 kel->cy = cy;
299 kel->cx = cx;
300 return 0;
301 }
302
303
304 /*!
305 * \brief kernelGetSum()
306 *
307 * \param[in] kel kernel
308 * \param[out] psum sum of all kernel values
309 * \return 0 if OK, 1 on error
310 */
311 l_ok
312 kernelGetSum(L_KERNEL *kel,
313 l_float32 *psum)
314 {
315 l_int32 sx, sy, i, j;
316
317 if (!psum)
318 return ERROR_INT("&sum not defined", __func__, 1);
319 *psum = 0.0;
320 if (!kel)
321 return ERROR_INT("kernel not defined", __func__, 1);
322
323 kernelGetParameters(kel, &sy, &sx, NULL, NULL);
324 for (i = 0; i < sy; i++) {
325 for (j = 0; j < sx; j++) {
326 *psum += kel->data[i][j];
327 }
328 }
329 return 0;
330 }
331
332
333 /*!
334 * \brief kernelGetMinMax()
335 *
336 * \param[in] kel kernel
337 * \param[out] pmin [optional] minimum value
338 * \param[out] pmax [optional] maximum value
339 * \return 0 if OK, 1 on error
340 */
341 l_ok
342 kernelGetMinMax(L_KERNEL *kel,
343 l_float32 *pmin,
344 l_float32 *pmax)
345 {
346 l_int32 sx, sy, i, j;
347 l_float32 val, minval, maxval;
348
349 if (!pmin && !pmax)
350 return ERROR_INT("neither &min nor &max defined", __func__, 1);
351 if (pmin) *pmin = 0.0;
352 if (pmax) *pmax = 0.0;
353 if (!kel)
354 return ERROR_INT("kernel not defined", __func__, 1);
355
356 kernelGetParameters(kel, &sy, &sx, NULL, NULL);
357 minval = 10000000.0;
358 maxval = -10000000.0;
359 for (i = 0; i < sy; i++) {
360 for (j = 0; j < sx; j++) {
361 val = kel->data[i][j];
362 if (val < minval)
363 minval = val;
364 if (val > maxval)
365 maxval = val;
366 }
367 }
368 if (pmin)
369 *pmin = minval;
370 if (pmax)
371 *pmax = maxval;
372
373 return 0;
374 }
375
376
377 /*----------------------------------------------------------------------*
378 * Normalize/Invert *
379 *----------------------------------------------------------------------*/
380 /*!
381 * \brief kernelNormalize()
382 *
383 * \param[in] kels source kel, to be normalized
384 * \param[in] normsum desired sum of elements in keld
385 * \return keld normalized version of kels, or NULL on error
386 * or if sum of elements is very close to 0)
387 *
388 * <pre>
389 * Notes:
390 * (1) If the sum of kernel elements is close to 0, do not
391 * try to calculate the normalized kernel. Instead,
392 * return a copy of the input kernel, with a warning.
393 * </pre>
394 */
395 L_KERNEL *
396 kernelNormalize(L_KERNEL *kels,
397 l_float32 normsum)
398 {
399 l_int32 i, j, sx, sy, cx, cy;
400 l_float32 sum, factor;
401 L_KERNEL *keld;
402
403 if (!kels)
404 return (L_KERNEL *)ERROR_PTR("kels not defined", __func__, NULL);
405
406 kernelGetSum(kels, &sum);
407 if (L_ABS(sum) < 0.00001) {
408 L_WARNING("null sum; not normalizing; returning a copy\n", __func__);
409 return kernelCopy(kels);
410 }
411
412 kernelGetParameters(kels, &sy, &sx, &cy, &cx);
413 if ((keld = kernelCreate(sy, sx)) == NULL)
414 return (L_KERNEL *)ERROR_PTR("keld not made", __func__, NULL);
415 keld->cy = cy;
416 keld->cx = cx;
417
418 factor = normsum / sum;
419 for (i = 0; i < sy; i++)
420 for (j = 0; j < sx; j++)
421 keld->data[i][j] = factor * kels->data[i][j];
422
423 return keld;
424 }
425
426
427 /*!
428 * \brief kernelInvert()
429 *
430 * \param[in] kels source kel, to be inverted
431 * \return keld spatially inverted, about the origin, or NULL on error
432 *
433 * <pre>
434 * Notes:
435 * (1) For convolution, the kernel is spatially inverted before
436 * a "correlation" operation is done between the kernel and the image.
437 * </pre>
438 */
439 L_KERNEL *
440 kernelInvert(L_KERNEL *kels)
441 {
442 l_int32 i, j, sx, sy, cx, cy;
443 L_KERNEL *keld;
444
445 if (!kels)
446 return (L_KERNEL *)ERROR_PTR("kels not defined", __func__, NULL);
447
448 kernelGetParameters(kels, &sy, &sx, &cy, &cx);
449 if ((keld = kernelCreate(sy, sx)) == NULL)
450 return (L_KERNEL *)ERROR_PTR("keld not made", __func__, NULL);
451 keld->cy = sy - 1 - cy;
452 keld->cx = sx - 1 - cx;
453
454 for (i = 0; i < sy; i++)
455 for (j = 0; j < sx; j++)
456 keld->data[i][j] = kels->data[sy - 1 - i][sx - 1 - j];
457
458 return keld;
459 }
460
461
462 /*----------------------------------------------------------------------*
463 * Helper function *
464 *----------------------------------------------------------------------*/
465 /*!
466 * \brief create2dFloatArray()
467 *
468 * \param[in] sy rows == height
469 * \param[in] sx columns == width
470 * \return doubly indexed array i.e., an array of sy row pointers,
471 * each of which points to an array of sx floats
472 *
473 * <pre>
474 * Notes:
475 * (1) The array[%sy][%sx] is indexed in standard "matrix notation",
476 * with the row index first.
477 * (2) The caller kernelCreate() limits the size to < 2^29 pixels.
478 * </pre>
479 */
480 l_float32 **
481 create2dFloatArray(l_int32 sy,
482 l_int32 sx)
483 {
484 l_int32 i;
485 l_float32 **array;
486
487 if (sx <= 0 || sx > MaxArraySize)
488 return (l_float32 **)ERROR_PTR("sx out of bounds", __func__, NULL);
489 if (sy <= 0 || sy > MaxArraySize)
490 return (l_float32 **)ERROR_PTR("sy out of bounds", __func__, NULL);
491
492 array = (l_float32 **)LEPT_CALLOC(sy, sizeof(l_float32 *));
493 for (i = 0; i < sy; i++)
494 array[i] = (l_float32 *)LEPT_CALLOC(sx, sizeof(l_float32));
495 return array;
496 }
497
498
499 /*----------------------------------------------------------------------*
500 * Kernel serialized I/O *
501 *----------------------------------------------------------------------*/
502 /*!
503 * \brief kernelRead()
504 *
505 * \param[in] fname filename
506 * \return kernel, or NULL on error
507 */
508 L_KERNEL *
509 kernelRead(const char *fname)
510 {
511 FILE *fp;
512 L_KERNEL *kel;
513
514 if (!fname)
515 return (L_KERNEL *)ERROR_PTR("fname not defined", __func__, NULL);
516
517 if ((fp = fopenReadStream(fname)) == NULL)
518 return (L_KERNEL *)ERROR_PTR_1("stream not opened",
519 fname, __func__, NULL);
520 if ((kel = kernelReadStream(fp)) == NULL) {
521 fclose(fp);
522 return (L_KERNEL *)ERROR_PTR_1("kel not returned",
523 fname, __func__, NULL);
524 }
525 fclose(fp);
526
527 return kel;
528 }
529
530
531 /*!
532 * \brief kernelReadStream()
533 *
534 * \param[in] fp file stream
535 * \return kernel, or NULL on error
536 */
537 L_KERNEL *
538 kernelReadStream(FILE *fp)
539 {
540 l_int32 sy, sx, cy, cx, i, j, ret, version, ignore;
541 L_KERNEL *kel;
542
543 if (!fp)
544 return (L_KERNEL *)ERROR_PTR("stream not defined", __func__, NULL);
545
546 ret = fscanf(fp, " Kernel Version %d\n", &version);
547 if (ret != 1)
548 return (L_KERNEL *)ERROR_PTR("not a kernel file", __func__, NULL);
549 if (version != KERNEL_VERSION_NUMBER)
550 return (L_KERNEL *)ERROR_PTR("invalid kernel version", __func__, NULL);
551
552 if (fscanf(fp, " sy = %d, sx = %d, cy = %d, cx = %d\n",
553 &sy, &sx, &cy, &cx) != 4)
554 return (L_KERNEL *)ERROR_PTR("dimensions not read", __func__, NULL);
555 if (sx > MaxArraySize || sy > MaxArraySize) {
556 L_ERROR("sx = %d or sy = %d > %d\n", __func__, sx, sy, MaxArraySize);
557 return NULL;
558 }
559 if ((kel = kernelCreate(sy, sx)) == NULL)
560 return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
561 kernelSetOrigin(kel, cy, cx);
562
563 for (i = 0; i < sy; i++) {
564 for (j = 0; j < sx; j++)
565 ignore = fscanf(fp, "%15f", &kel->data[i][j]);
566 ignore = fscanf(fp, "\n");
567 }
568 ignore = fscanf(fp, "\n");
569
570 return kel;
571 }
572
573
574 /*!
575 * \brief kernelWrite()
576 *
577 * \param[in] fname output file
578 * \param[in] kel kernel
579 * \return 0 if OK, 1 on error
580 */
581 l_ok
582 kernelWrite(const char *fname,
583 L_KERNEL *kel)
584 {
585 FILE *fp;
586
587 if (!fname)
588 return ERROR_INT("fname not defined", __func__, 1);
589 if (!kel)
590 return ERROR_INT("kel not defined", __func__, 1);
591
592 if ((fp = fopenWriteStream(fname, "wb")) == NULL)
593 return ERROR_INT_1("stream not opened", fname, __func__, 1);
594 kernelWriteStream(fp, kel);
595 fclose(fp);
596
597 return 0;
598 }
599
600
601 /*!
602 * \brief kernelWriteStream()
603 *
604 * \param[in] fp file stream
605 * \param[in] kel
606 * \return 0 if OK, 1 on error
607 */
608 l_ok
609 kernelWriteStream(FILE *fp,
610 L_KERNEL *kel)
611 {
612 l_int32 sx, sy, cx, cy, i, j;
613
614 if (!fp)
615 return ERROR_INT("stream not defined", __func__, 1);
616 if (!kel)
617 return ERROR_INT("kel not defined", __func__, 1);
618 kernelGetParameters(kel, &sy, &sx, &cy, &cx);
619
620 fprintf(fp, " Kernel Version %d\n", KERNEL_VERSION_NUMBER);
621 fprintf(fp, " sy = %d, sx = %d, cy = %d, cx = %d\n", sy, sx, cy, cx);
622 for (i = 0; i < sy; i++) {
623 for (j = 0; j < sx; j++)
624 fprintf(fp, "%15.4f", kel->data[i][j]);
625 fprintf(fp, "\n");
626 }
627 fprintf(fp, "\n");
628
629 return 0;
630 }
631
632
633 /*----------------------------------------------------------------------*
634 * Making a kernel from a compiled string *
635 *----------------------------------------------------------------------*/
636 /*!
637 * \brief kernelCreateFromString()
638 *
639 * \param[in] h, w height, width
640 * \param[in] cy, cx origin
641 * \param[in] kdata
642 * \return kernel of the given size, or NULL on error
643 *
644 * <pre>
645 * Notes:
646 * (1) The data is an array of chars, in row-major order, giving
647 * space separated integers in the range [-255 ... 255].
648 * (2) The only other formatting limitation is that you must
649 * leave space between the last number in each row and
650 * the double-quote. If possible, it's also nice to have each
651 * line in the string represent a line in the kernel; e.g.,
652 * static const char *kdata =
653 * " 20 50 20 "
654 * " 70 140 70 "
655 * " 20 50 20 ";
656 * </pre>
657 */
658 L_KERNEL *
659 kernelCreateFromString(l_int32 h,
660 l_int32 w,
661 l_int32 cy,
662 l_int32 cx,
663 const char *kdata)
664 {
665 l_int32 n, i, j, index;
666 l_float32 val;
667 L_KERNEL *kel;
668 NUMA *na;
669
670 if (h < 1)
671 return (L_KERNEL *)ERROR_PTR("height must be > 0", __func__, NULL);
672 if (w < 1)
673 return (L_KERNEL *)ERROR_PTR("width must be > 0", __func__, NULL);
674 if (cy < 0 || cy >= h)
675 return (L_KERNEL *)ERROR_PTR("cy invalid", __func__, NULL);
676 if (cx < 0 || cx >= w)
677 return (L_KERNEL *)ERROR_PTR("cx invalid", __func__, NULL);
678
679 kel = kernelCreate(h, w);
680 kernelSetOrigin(kel, cy, cx);
681 na = parseStringForNumbers(kdata, " \t\n");
682 n = numaGetCount(na);
683 if (n != w * h) {
684 kernelDestroy(&kel);
685 numaDestroy(&na);
686 lept_stderr("w = %d, h = %d, num ints = %d\n", w, h, n);
687 return (L_KERNEL *)ERROR_PTR("invalid integer data", __func__, NULL);
688 }
689
690 index = 0;
691 for (i = 0; i < h; i++) {
692 for (j = 0; j < w; j++) {
693 numaGetFValue(na, index, &val);
694 kernelSetElement(kel, i, j, val);
695 index++;
696 }
697 }
698
699 numaDestroy(&na);
700 return kel;
701 }
702
703
704 /*----------------------------------------------------------------------*
705 * Making a kernel from a simple file format *
706 *----------------------------------------------------------------------*/
707 /*!
708 * \brief kernelCreateFromFile()
709 *
710 * \param[in] filename
711 * \return kernel, or NULL on error
712 *
713 * <pre>
714 * Notes:
715 * (1) The file contains, in the following order:
716 * ~ Any number of comment lines starting with '#' are ignored
717 * ~ The height and width of the kernel
718 * ~ The y and x values of the kernel origin
719 * ~ The kernel data, formatted as lines of numbers (integers
720 * or floats) for the kernel values in row-major order,
721 * and with no other punctuation.
722 * (Note: this differs from kernelCreateFromString(),
723 * where each line must begin and end with a double-quote
724 * to tell the compiler it's part of a string.)
725 * ~ The kernel specification ends when a blank line,
726 * a comment line, or the end of file is reached.
727 * (2) All lines must be left-justified.
728 * (3) See kernelCreateFromString() for a description of the string
729 * format for the kernel data. As an example, here are the lines
730 * of a valid kernel description file In the file, all lines
731 * are left-justified:
732 * \code
733 * # small 3x3 kernel
734 * 3 3
735 * 1 1
736 * 25.5 51 24.3
737 * 70.2 146.3 73.4
738 * 20 50.9 18.4
739 * \endcode
740 * </pre>
741 */
742 L_KERNEL *
743 kernelCreateFromFile(const char *filename)
744 {
745 char *filestr, *line;
746 l_int32 nlines, i, j, first, index, w, h, cx, cy, n;
747 l_float32 val;
748 size_t size;
749 NUMA *na, *nat;
750 SARRAY *sa;
751 L_KERNEL *kel;
752
753 if (!filename)
754 return (L_KERNEL *)ERROR_PTR("filename not defined", __func__, NULL);
755
756 if ((filestr = (char *)l_binaryRead(filename, &size)) == NULL)
757 return (L_KERNEL *)ERROR_PTR_1("file not found",
758 filename, __func__, NULL);
759 if (size == 0) {
760 LEPT_FREE(filestr);
761 return (L_KERNEL *)ERROR_PTR_1("file is empty",
762 filename, __func__, NULL);
763 }
764
765 sa = sarrayCreateLinesFromString(filestr, 1);
766 LEPT_FREE(filestr);
767 nlines = sarrayGetCount(sa);
768
769 /* Find the first data line. */
770 for (i = 0, first = 0; i < nlines; i++) {
771 line = sarrayGetString(sa, i, L_NOCOPY);
772 if (line[0] != '#') {
773 first = i;
774 break;
775 }
776 }
777
778 /* Find the kernel dimensions and origin location. */
779 line = sarrayGetString(sa, first, L_NOCOPY);
780 if (sscanf(line, "%d %d", &h, &w) != 2) {
781 sarrayDestroy(&sa);
782 return (L_KERNEL *)ERROR_PTR("error reading h,w", __func__, NULL);
783 }
784 if (h > MaxArraySize || w > MaxArraySize) {
785 L_ERROR("h = %d or w = %d > %d\n", __func__, h, w, MaxArraySize);
786 sarrayDestroy(&sa);
787 return NULL;
788 }
789 line = sarrayGetString(sa, first + 1, L_NOCOPY);
790 if (sscanf(line, "%d %d", &cy, &cx) != 2) {
791 sarrayDestroy(&sa);
792 return (L_KERNEL *)ERROR_PTR("error reading cy,cx", __func__, NULL);
793 }
794
795 /* Extract the data. This ends when we reach eof, or when we
796 * encounter a line of data that is either a null string or
797 * contains just a newline. */
798 na = numaCreate(0);
799 for (i = first + 2; i < nlines; i++) {
800 line = sarrayGetString(sa, i, L_NOCOPY);
801 if (line[0] == '\0' || line[0] == '\n' || line[0] == '#')
802 break;
803 nat = parseStringForNumbers(line, " \t\n");
804 numaJoin(na, nat, 0, -1);
805 numaDestroy(&nat);
806 }
807 sarrayDestroy(&sa);
808
809 n = numaGetCount(na);
810 if (n != w * h) {
811 numaDestroy(&na);
812 lept_stderr("w = %d, h = %d, num ints = %d\n", w, h, n);
813 return (L_KERNEL *)ERROR_PTR("invalid integer data", __func__, NULL);
814 }
815
816 kel = kernelCreate(h, w);
817 kernelSetOrigin(kel, cy, cx);
818 index = 0;
819 for (i = 0; i < h; i++) {
820 for (j = 0; j < w; j++) {
821 numaGetFValue(na, index, &val);
822 kernelSetElement(kel, i, j, val);
823 index++;
824 }
825 }
826
827 numaDestroy(&na);
828 return kel;
829 }
830
831
832 /*----------------------------------------------------------------------*
833 * Making a kernel from a Pix *
834 *----------------------------------------------------------------------*/
835 /*!
836 * \brief kernelCreateFromPix()
837 *
838 * \param[in] pix
839 * \param[in] cy, cx origin of kernel
840 * \return kernel, or NULL on error
841 *
842 * <pre>
843 * Notes:
844 * (1) The origin must be positive and within the dimensions of the pix.
845 * </pre>
846 */
847 L_KERNEL *
848 kernelCreateFromPix(PIX *pix,
849 l_int32 cy,
850 l_int32 cx)
851 {
852 l_int32 i, j, w, h, d;
853 l_uint32 val;
854 L_KERNEL *kel;
855
856 if (!pix)
857 return (L_KERNEL *)ERROR_PTR("pix not defined", __func__, NULL);
858 pixGetDimensions(pix, &w, &h, &d);
859 if (d != 8)
860 return (L_KERNEL *)ERROR_PTR("pix not 8 bpp", __func__, NULL);
861 if (cy < 0 || cx < 0 || cy >= h || cx >= w)
862 return (L_KERNEL *)ERROR_PTR("(cy, cx) invalid", __func__, NULL);
863
864 kel = kernelCreate(h, w);
865 kernelSetOrigin(kel, cy, cx);
866 for (i = 0; i < h; i++) {
867 for (j = 0; j < w; j++) {
868 pixGetPixel(pix, j, i, &val);
869 kernelSetElement(kel, i, j, (l_float32)val);
870 }
871 }
872
873 return kel;
874 }
875
876
877 /*----------------------------------------------------------------------*
878 * Display a kernel in a pix *
879 *----------------------------------------------------------------------*/
880 /*!
881 * \brief kernelDisplayInPix()
882 *
883 * \param[in] kel kernel
884 * \param[in] size of grid interiors; odd; either 1 or a minimum size
885 * of 17 is enforced
886 * \param[in] gthick grid thickness; either 0 or a minimum size of 2
887 * is enforced
888 * \return pix display of kernel, or NULL on error
889 *
890 * <pre>
891 * Notes:
892 * (1) This gives a visual representation of a kernel.
893 * (2) There are two modes of display:
894 * (a) Grid lines of minimum width 2, surrounding regions
895 * representing kernel elements of minimum size 17,
896 * with a "plus" mark at the kernel origin, or
897 * (b) A pix without grid lines and using 1 pixel per kernel element.
898 * (3) For both cases, the kernel absolute value is displayed,
899 * normalized such that the maximum absolute value is 255.
900 * (4) Large 2D separable kernels should be used for convolution
901 * with two 1D kernels. However, for the bilateral filter,
902 * the computation time is independent of the size of the
903 * 2D content kernel.
904 * </pre>
905 */
906 PIX *
907 kernelDisplayInPix(L_KERNEL *kel,
908 l_int32 size,
909 l_int32 gthick)
910 {
911 l_int32 i, j, w, h, sx, sy, cx, cy, width, x0, y0;
912 l_int32 normval;
913 l_float32 minval, maxval, max, val, norm;
914 PIX *pixd, *pixt0, *pixt1;
915
916 if (!kel)
917 return (PIX *)ERROR_PTR("kernel not defined", __func__, NULL);
918
919 /* Normalize the max value to be 255 for display */
920 kernelGetParameters(kel, &sy, &sx, &cy, &cx);
921 kernelGetMinMax(kel, &minval, &maxval);
922 max = L_MAX(maxval, -minval);
923 if (max == 0.0)
924 return (PIX *)ERROR_PTR("kernel elements all 0.0", __func__, NULL);
925 norm = 255.f / (l_float32)max;
926
927 /* Handle the 1 element/pixel case; typically with large kernels */
928 if (size == 1 && gthick == 0) {
929 pixd = pixCreate(sx, sy, 8);
930 for (i = 0; i < sy; i++) {
931 for (j = 0; j < sx; j++) {
932 kernelGetElement(kel, i, j, &val);
933 normval = (l_int32)(norm * L_ABS(val));
934 pixSetPixel(pixd, j, i, normval);
935 }
936 }
937 return pixd;
938 }
939
940 /* Enforce the constraints for the grid line version */
941 if (size < 17) {
942 L_WARNING("size < 17; setting to 17\n", __func__);
943 size = 17;
944 }
945 if (size % 2 == 0)
946 size++;
947 if (gthick < 2) {
948 L_WARNING("grid thickness < 2; setting to 2\n", __func__);
949 gthick = 2;
950 }
951
952 w = size * sx + gthick * (sx + 1);
953 h = size * sy + gthick * (sy + 1);
954 pixd = pixCreate(w, h, 8);
955
956 /* Generate grid lines */
957 for (i = 0; i <= sy; i++)
958 pixRenderLine(pixd, 0, gthick / 2 + i * (size + gthick),
959 w - 1, gthick / 2 + i * (size + gthick),
960 gthick, L_SET_PIXELS);
961 for (j = 0; j <= sx; j++)
962 pixRenderLine(pixd, gthick / 2 + j * (size + gthick), 0,
963 gthick / 2 + j * (size + gthick), h - 1,
964 gthick, L_SET_PIXELS);
965
966 /* Generate mask for each element */
967 pixt0 = pixCreate(size, size, 1);
968 pixSetAll(pixt0);
969
970 /* Generate crossed lines for origin pattern */
971 pixt1 = pixCreate(size, size, 1);
972 width = size / 8;
973 pixRenderLine(pixt1, size / 2, (l_int32)(0.12 * size),
974 size / 2, (l_int32)(0.88 * size),
975 width, L_SET_PIXELS);
976 pixRenderLine(pixt1, (l_int32)(0.15 * size), size / 2,
977 (l_int32)(0.85 * size), size / 2,
978 width, L_FLIP_PIXELS);
979 pixRasterop(pixt1, size / 2 - width, size / 2 - width,
980 2 * width, 2 * width, PIX_NOT(PIX_DST), NULL, 0, 0);
981
982 /* Paste the patterns in */
983 y0 = gthick;
984 for (i = 0; i < sy; i++) {
985 x0 = gthick;
986 for (j = 0; j < sx; j++) {
987 kernelGetElement(kel, i, j, &val);
988 normval = (l_int32)(norm * L_ABS(val));
989 pixSetMaskedGeneral(pixd, pixt0, normval, x0, y0);
990 if (i == cy && j == cx)
991 pixPaintThroughMask(pixd, pixt1, x0, y0, 255 - normval);
992 x0 += size + gthick;
993 }
994 y0 += size + gthick;
995 }
996
997 pixDestroy(&pixt0);
998 pixDestroy(&pixt1);
999 return pixd;
1000 }
1001
1002
1003 /*------------------------------------------------------------------------*
1004 * Parse string to extract numbers *
1005 *------------------------------------------------------------------------*/
1006 /*!
1007 * \brief parseStringForNumbers()
1008 *
1009 * \param[in] str string containing numbers; not changed
1010 * \param[in] seps string of characters that can be used between ints
1011 * \return numa of numbers found, or NULL on error
1012 *
1013 * <pre>
1014 * Notes:
1015 * (1) The numbers can be ints or floats.
1016 * </pre>
1017 */
1018 NUMA *
1019 parseStringForNumbers(const char *str,
1020 const char *seps)
1021 {
1022 char *newstr, *head;
1023 char *tail = NULL;
1024 l_float32 val;
1025 NUMA *na;
1026
1027 if (!str)
1028 return (NUMA *)ERROR_PTR("str not defined", __func__, NULL);
1029
1030 newstr = stringNew(str); /* to enforce const-ness of str */
1031 na = numaCreate(0);
1032 head = strtokSafe(newstr, seps, &tail);
1033 val = atof(head);
1034 numaAddNumber(na, val);
1035 LEPT_FREE(head);
1036 while ((head = strtokSafe(NULL, seps, &tail)) != NULL) {
1037 val = atof(head);
1038 numaAddNumber(na, val);
1039 LEPT_FREE(head);
1040 }
1041
1042 LEPT_FREE(newstr);
1043 return na;
1044 }
1045
1046
1047 /*------------------------------------------------------------------------*
1048 * Simple parametric kernels *
1049 *------------------------------------------------------------------------*/
1050 /*!
1051 * \brief makeFlatKernel()
1052 *
1053 * \param[in] height, width
1054 * \param[in] cy, cx origin of kernel
1055 * \return kernel, or NULL on error
1056 *
1057 * <pre>
1058 * Notes:
1059 * (1) This is the same low-pass filtering kernel that is used
1060 * in the block convolution functions.
1061 * (2) The kernel origin (%cy, %cx) is typically placed as near
1062 * the center of the kernel as possible. If height and
1063 * width are odd, then using %cy = height / 2 and
1064 * %cx = width / 2 places the origin at the exact center.
1065 * (3) This returns a normalized kernel.
1066 * </pre>
1067 */
1068 L_KERNEL *
1069 makeFlatKernel(l_int32 height,
1070 l_int32 width,
1071 l_int32 cy,
1072 l_int32 cx)
1073 {
1074 l_int32 i, j;
1075 l_float32 normval;
1076 L_KERNEL *kel;
1077
1078 if ((kel = kernelCreate(height, width)) == NULL)
1079 return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
1080 kernelSetOrigin(kel, cy, cx);
1081 normval = 1.0f / (l_float32)(height * width);
1082 for (i = 0; i < height; i++) {
1083 for (j = 0; j < width; j++) {
1084 kernelSetElement(kel, i, j, normval);
1085 }
1086 }
1087
1088 return kel;
1089 }
1090
1091
1092 /*!
1093 * \brief makeGaussianKernel()
1094 *
1095 * \param[in] halfh sy = 2 * halfh + 1
1096 * \param[in] halfw sx = 2 * halfw + 1
1097 * \param[in] stdev standard deviation
1098 * \param[in] max value at (cx,cy)
1099 * \return kernel, or NULL on error
1100 *
1101 * <pre>
1102 * Notes:
1103 * (1) The kernel size (sx, sy) = (2 * %halfw + 1, 2 * %halfh + 1)
1104 * (2) The kernel center (cx, cy) = (%halfw, %halfh).
1105 * (3) %halfw and %halfh are typically equal, and
1106 * are typically several times larger than the standard deviation.
1107 * (4) If pixConvolve() is invoked with normalization (the sum of
1108 * kernel elements = 1.0), use 1.0 for max (or any number that's
1109 * not too small or too large).
1110 * </pre>
1111 */
1112 L_KERNEL *
1113 makeGaussianKernel(l_int32 halfh,
1114 l_int32 halfw,
1115 l_float32 stdev,
1116 l_float32 max)
1117 {
1118 l_int32 sx, sy, i, j;
1119 l_float32 val;
1120 L_KERNEL *kel;
1121
1122 sx = 2 * halfw + 1;
1123 sy = 2 * halfh + 1;
1124 if ((kel = kernelCreate(sy, sx)) == NULL)
1125 return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
1126 kernelSetOrigin(kel, halfh, halfw);
1127 for (i = 0; i < sy; i++) {
1128 for (j = 0; j < sx; j++) {
1129 val = expf(-(l_float32)((i - halfh) * (i - halfh) +
1130 (j - halfw) * (j - halfw)) /
1131 (2. * stdev * stdev));
1132 kernelSetElement(kel, i, j, max * val);
1133 }
1134 }
1135
1136 return kel;
1137 }
1138
1139
1140 /*!
1141 * \brief makeGaussianKernelSep()
1142 *
1143 * \param[in] halfh sy = 2 * halfh + 1
1144 * \param[in] halfw sx = 2 * halfw + 1
1145 * \param[in] stdev standard deviation
1146 * \param[in] max value at (cx,cy)
1147 * \param[out] pkelx x part of kernel
1148 * \param[out] pkely y part of kernel
1149 * \return 0 if OK, 1 on error
1150 *
1151 * <pre>
1152 * Notes:
1153 * (1) See makeGaussianKernel() for description of input parameters.
1154 * (2) These kernels are constructed so that the result of both
1155 * normalized and un-normalized convolution will be the same
1156 * as when convolving with pixConvolve() using the full kernel.
1157 * (3) The trick for the un-normalized convolution is to have the
1158 * product of the two kernel elements at (cx,cy) be equal to %max,
1159 * not max**2. That's why %max for kely is 1.0. If instead
1160 * we use sqrt(%max) for both, the results are slightly less
1161 * accurate, when compared to using the full kernel in
1162 * makeGaussianKernel().
1163 * </pre>
1164 */
1165 l_ok
1166 makeGaussianKernelSep(l_int32 halfh,
1167 l_int32 halfw,
1168 l_float32 stdev,
1169 l_float32 max,
1170 L_KERNEL **pkelx,
1171 L_KERNEL **pkely)
1172 {
1173 if (!pkelx || !pkely)
1174 return ERROR_INT("&kelx and &kely not defined", __func__, 1);
1175
1176 *pkelx = makeGaussianKernel(0, halfw, stdev, max);
1177 *pkely = makeGaussianKernel(halfh, 0, stdev, 1.0);
1178 return 0;
1179 }
1180
1181
1182 /*!
1183 * \brief makeDoGKernel()
1184 *
1185 * \param[in] halfh sy = 2 * halfh + 1
1186 * \param[in] halfw sx = 2 * halfw + 1
1187 * \param[in] stdev standard deviation of narrower gaussian
1188 * \param[in] ratio of stdev for wide filter to stdev for narrow one
1189 * \return kernel, or NULL on error
1190 *
1191 * <pre>
1192 * Notes:
1193 * (1) The DoG (difference of gaussians) is a wavelet mother
1194 * function with null total sum. By subtracting two blurred
1195 * versions of the image, it acts as a bandpass filter for
1196 * frequencies passed by the narrow gaussian but stopped
1197 * by the wide one.See:
1198 * http://en.wikipedia.org/wiki/Difference_of_Gaussians
1199 * (2) The kernel size (sx, sy) = (2 * halfw + 1, 2 * halfh + 1).
1200 * (3) The kernel center (cx, cy) = (halfw, halfh).
1201 * (4) %halfw and %halfh are typically equal, and are typically
1202 * several times larger than the standard deviation.
1203 * (5) %ratio is the ratio of standard deviations of the wide
1204 * to narrow gaussian. It must be >= 1.0; 1.0 is a no-op.
1205 * (6) Because the kernel is a null sum, it must be invoked without
1206 * normalization in pixConvolve().
1207 * </pre>
1208 */
1209 L_KERNEL *
1210 makeDoGKernel(l_int32 halfh,
1211 l_int32 halfw,
1212 l_float32 stdev,
1213 l_float32 ratio)
1214 {
1215 l_int32 sx, sy, i, j;
1216 l_float32 pi, squaredist, highnorm, lownorm, val;
1217 L_KERNEL *kel;
1218
1219 sx = 2 * halfw + 1;
1220 sy = 2 * halfh + 1;
1221 if ((kel = kernelCreate(sy, sx)) == NULL)
1222 return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
1223 kernelSetOrigin(kel, halfh, halfw);
1224
1225 pi = 3.1415926535f;
1226 for (i = 0; i < sy; i++) {
1227 for (j = 0; j < sx; j++) {
1228 squaredist = (l_float32)((i - halfh) * (i - halfh) +
1229 (j - halfw) * (j - halfw));
1230 highnorm = 1.f / (2 * stdev * stdev);
1231 lownorm = highnorm / (ratio * ratio);
1232 val = (highnorm / pi) * expf(-(highnorm * squaredist))
1233 - (lownorm / pi) * expf(-(lownorm * squaredist));
1234 kernelSetElement(kel, i, j, val);
1235 }
1236 }
1237
1238 return kel;
1239 }