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