Mercurial > hgrepos > Python2 > PyMuPDF

comparison mupdf-source/thirdparty/leptonica/src/ptafunc1.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 * \file ptafunc1.c
29 * <pre>
30 *
31 * --------------------------------------
32 * This file has these Pta utilities:
33 * - simple rearrangements
34 * - geometric analysis
35 * - min/max and filtering
36 * - least squares fitting
37 * - interconversions with Pix and Numa
38 * - display into a pix
39 * --------------------------------------
40 *
41 * Simple rearrangements
42 * PTA *ptaSubsample()
43 * l_int32 ptaJoin()
44 * l_int32 ptaaJoin()
45 * PTA *ptaReverse()
46 * PTA *ptaTranspose()
47 * PTA *ptaCyclicPerm()
48 * PTA *ptaSelectRange()
49 *
50 * Geometric
51 * BOX *ptaGetBoundingRegion()
52 * l_int32 *ptaGetRange()
53 * PTA *ptaGetInsideBox()
54 * PTA *pixFindCornerPixels()
55 * l_int32 ptaContainsPt()
56 * l_int32 ptaTestIntersection()
57 * PTA *ptaTransform()
58 * l_int32 ptaPtInsidePolygon()
59 * l_float32 l_angleBetweenVectors()
60 * l_int32 ptaPolygonIsConvex()
61 *
62 * Min/max and filtering
63 * l_int32 ptaGetMinMax()
64 * PTA *ptaSelectByValue()
65 * PTA *ptaCropToMask()
66 *
67 * Least Squares Fit
68 * l_int32 ptaGetLinearLSF()
69 * l_int32 ptaGetQuadraticLSF()
70 * l_int32 ptaGetCubicLSF()
71 * l_int32 ptaGetQuarticLSF()
72 * l_int32 ptaNoisyLinearLSF()
73 * l_int32 ptaNoisyQuadraticLSF()
74 * l_int32 applyLinearFit()
75 * l_int32 applyQuadraticFit()
76 * l_int32 applyCubicFit()
77 * l_int32 applyQuarticFit()
78 *
79 * Interconversions with Pix
80 * l_int32 pixPlotAlongPta()
81 * PTA *ptaGetPixelsFromPix()
82 * PIX *pixGenerateFromPta()
83 * PTA *ptaGetBoundaryPixels()
84 * PTAA *ptaaGetBoundaryPixels()
85 * PTAA *ptaaIndexLabeledPixels()
86 * PTA *ptaGetNeighborPixLocs()
87 *
88 * Interconversion with Numa
89 * PTA *numaConvertToPta1()
90 * PTA *numaConvertToPta2()
91 * l_int32 ptaConvertToNuma()
92 *
93 * Display Pta and Ptaa
94 * PIX *pixDisplayPta()
95 * PIX *pixDisplayPtaaPattern()
96 * PIX *pixDisplayPtaPattern()
97 * PTA *ptaReplicatePattern()
98 * PIX *pixDisplayPtaa()
99 * </pre>
100 */
101
102 #ifdef HAVE_CONFIG_H
103 #include <config_auto.h>
104 #endif /* HAVE_CONFIG_H */
105
106 #include <math.h>
107 #include "allheaders.h"
108 #include "pix_internal.h"
109
110 #ifndef M_PI
111 #define M_PI 3.14159265358979323846
112 #endif /* M_PI */
113
114 /*---------------------------------------------------------------------*
115 * Simple rearrangements *
116 *---------------------------------------------------------------------*/
117 /*!
118 * \brief ptaSubsample()
119 *
120 * \param[in] ptas
121 * \param[in] subfactor subsample factor, >= 1
122 * \return ptad evenly sampled pt values from ptas, or NULL on error
123 */
124 PTA *
125 ptaSubsample(PTA *ptas,
126 l_int32 subfactor)
127 {
128 l_int32 n, i;
129 l_float32 x, y;
130 PTA *ptad;
131
132 if (!ptas)
133 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
134 if (subfactor < 1)
135 return (PTA *)ERROR_PTR("subfactor < 1", __func__, NULL);
136
137 ptad = ptaCreate(0);
138 n = ptaGetCount(ptas);
139 for (i = 0; i < n; i++) {
140 if (i % subfactor != 0) continue;
141 ptaGetPt(ptas, i, &x, &y);
142 ptaAddPt(ptad, x, y);
143 }
144
145 return ptad;
146 }
147
148
149 /*!
150 * \brief ptaJoin()
151 *
152 * \param[in] ptad dest pta; add to this one
153 * \param[in] ptas source pta; add from this one
154 * \param[in] istart starting index in ptas
155 * \param[in] iend ending index in ptas; use -1 to cat all
156 * \return 0 if OK, 1 on error
157 *
158 * <pre>
159 * Notes:
160 * (1) istart < 0 is taken to mean 'read from the start' (istart = 0)
161 * (2) iend < 0 means 'read to the end'
162 * (3) if ptas == NULL, this is a no-op
163 * </pre>
164 */
165 l_ok
166 ptaJoin(PTA *ptad,
167 PTA *ptas,
168 l_int32 istart,
169 l_int32 iend)
170 {
171 l_int32 n, i, x, y;
172
173 if (!ptad)
174 return ERROR_INT("ptad not defined", __func__, 1);
175 if (!ptas)
176 return 0;
177
178 if (istart < 0)
179 istart = 0;
180 n = ptaGetCount(ptas);
181 if (iend < 0 || iend >= n)
182 iend = n - 1;
183 if (istart > iend)
184 return ERROR_INT("istart > iend; no pts", __func__, 1);
185
186 for (i = istart; i <= iend; i++) {
187 ptaGetIPt(ptas, i, &x, &y);
188 if (ptaAddPt(ptad, x, y) == 1) {
189 L_ERROR("failed to add pt at i = %d\n", __func__, i);
190 return 1;
191 }
192 }
193 return 0;
194 }
195
196
197 /*!
198 * \brief ptaaJoin()
199 *
200 * \param[in] ptaad dest ptaa; add to this one
201 * \param[in] ptaas source ptaa; add from this one
202 * \param[in] istart starting index in ptaas
203 * \param[in] iend ending index in ptaas; use -1 to cat all
204 * \return 0 if OK, 1 on error
205 *
206 * <pre>
207 * Notes:
208 * (1) istart < 0 is taken to mean 'read from the start' (istart = 0)
209 * (2) iend < 0 means 'read to the end'
210 * (3) if ptas == NULL, this is a no-op
211 * </pre>
212 */
213 l_ok
214 ptaaJoin(PTAA *ptaad,
215 PTAA *ptaas,
216 l_int32 istart,
217 l_int32 iend)
218 {
219 l_int32 n, i;
220 PTA *pta;
221
222 if (!ptaad)
223 return ERROR_INT("ptaad not defined", __func__, 1);
224 if (!ptaas)
225 return 0;
226
227 if (istart < 0)
228 istart = 0;
229 n = ptaaGetCount(ptaas);
230 if (iend < 0 || iend >= n)
231 iend = n - 1;
232 if (istart > iend)
233 return ERROR_INT("istart > iend; no pts", __func__, 1);
234
235 for (i = istart; i <= iend; i++) {
236 pta = ptaaGetPta(ptaas, i, L_CLONE);
237 ptaaAddPta(ptaad, pta, L_INSERT);
238 }
239
240 return 0;
241 }
242
243
244 /*!
245 * \brief ptaReverse()
246 *
247 * \param[in] ptas
248 * \param[in] type 0 for float values; 1 for integer values
249 * \return ptad reversed pta, or NULL on error
250 */
251 PTA *
252 ptaReverse(PTA *ptas,
253 l_int32 type)
254 {
255 l_int32 n, i, ix, iy;
256 l_float32 x, y;
257 PTA *ptad;
258
259 if (!ptas)
260 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
261
262 n = ptaGetCount(ptas);
263 if ((ptad = ptaCreate(n)) == NULL)
264 return (PTA *)ERROR_PTR("ptad not made", __func__, NULL);
265 for (i = n - 1; i >= 0; i--) {
266 if (type == 0) {
267 ptaGetPt(ptas, i, &x, &y);
268 ptaAddPt(ptad, x, y);
269 } else { /* type == 1 */
270 ptaGetIPt(ptas, i, &ix, &iy);
271 ptaAddPt(ptad, ix, iy);
272 }
273 }
274
275 return ptad;
276 }
277
278
279 /*!
280 * \brief ptaTranspose()
281 *
282 * \param[in] ptas
283 * \return ptad with x and y values swapped, or NULL on error
284 */
285 PTA *
286 ptaTranspose(PTA *ptas)
287 {
288 l_int32 n, i;
289 l_float32 x, y;
290 PTA *ptad;
291
292 if (!ptas)
293 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
294
295 n = ptaGetCount(ptas);
296 if ((ptad = ptaCreate(n)) == NULL)
297 return (PTA *)ERROR_PTR("ptad not made", __func__, NULL);
298 for (i = 0; i < n; i++) {
299 ptaGetPt(ptas, i, &x, &y);
300 ptaAddPt(ptad, y, x);
301 }
302
303 return ptad;
304 }
305
306
307 /*!
308 * \brief ptaCyclicPerm()
309 *
310 * \param[in] ptas
311 * \param[in] xs, ys start point; must be in ptas
312 * \return ptad cyclic permutation, starting and ending at (xs, ys,
313 * or NULL on error
314 *
315 * <pre>
316 * Notes:
317 * (1) Check to insure that (a) ptas is a closed path where
318 * the first and last points are identical, and (b) the
319 * resulting pta also starts and ends on the same point
320 * (which in this case is (xs, ys).
321 * </pre>
322 */
323 PTA *
324 ptaCyclicPerm(PTA *ptas,
325 l_int32 xs,
326 l_int32 ys)
327 {
328 l_int32 n, i, x, y, j, index, state;
329 l_int32 x1, y1, x2, y2;
330 PTA *ptad;
331
332 if (!ptas)
333 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
334
335 n = ptaGetCount(ptas);
336
337 /* Verify input data */
338 ptaGetIPt(ptas, 0, &x1, &y1);
339 ptaGetIPt(ptas, n - 1, &x2, &y2);
340 if (x1 != x2 || y1 != y2)
341 return (PTA *)ERROR_PTR("start and end pts not same", __func__, NULL);
342 state = L_NOT_FOUND;
343 for (i = 0; i < n; i++) {
344 ptaGetIPt(ptas, i, &x, &y);
345 if (x == xs && y == ys) {
346 state = L_FOUND;
347 break;
348 }
349 }
350 if (state == L_NOT_FOUND)
351 return (PTA *)ERROR_PTR("start pt not in ptas", __func__, NULL);
352
353 if ((ptad = ptaCreate(n)) == NULL)
354 return (PTA *)ERROR_PTR("ptad not made", __func__, NULL);
355 for (j = 0; j < n - 1; j++) {
356 if (i + j < n - 1)
357 index = i + j;
358 else
359 index = (i + j + 1) % n;
360 ptaGetIPt(ptas, index, &x, &y);
361 ptaAddPt(ptad, x, y);
362 }
363 ptaAddPt(ptad, xs, ys);
364
365 return ptad;
366 }
367
368
369 /*!
370 * \brief ptaSelectRange()
371 *
372 * \param[in] ptas
373 * \param[in] first use 0 to select from the beginning
374 * \param[in] last use -1 to select to the end
375 * \return ptad, or NULL on error
376 */
377 PTA *
378 ptaSelectRange(PTA *ptas,
379 l_int32 first,
380 l_int32 last)
381 {
382 l_int32 n, npt, i;
383 l_float32 x, y;
384 PTA *ptad;
385
386 if (!ptas)
387 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
388 if ((n = ptaGetCount(ptas)) == 0) {
389 L_WARNING("ptas is empty\n", __func__);
390 return ptaCopy(ptas);
391 }
392 first = L_MAX(0, first);
393 if (last < 0) last = n - 1;
394 if (first >= n)
395 return (PTA *)ERROR_PTR("invalid first", __func__, NULL);
396 if (last >= n) {
397 L_WARNING("last = %d is beyond max index = %d; adjusting\n",
398 __func__, last, n - 1);
399 last = n - 1;
400 }
401 if (first > last)
402 return (PTA *)ERROR_PTR("first > last", __func__, NULL);
403
404 npt = last - first + 1;
405 ptad = ptaCreate(npt);
406 for (i = first; i <= last; i++) {
407 ptaGetPt(ptas, i, &x, &y);
408 ptaAddPt(ptad, x, y);
409 }
410 return ptad;
411 }
412
413
414 /*---------------------------------------------------------------------*
415 * Geometric *
416 *---------------------------------------------------------------------*/
417 /*!
418 * \brief ptaGetBoundingRegion()
419 *
420 * \param[in] pta
421 * \return box, or NULL on error
422 *
423 * <pre>
424 * Notes:
425 * (1) This is used when the pta represents a set of points in
426 * a two-dimensional image. It returns the box of minimum
427 * size containing the pts in the pta.
428 * </pre>
429 */
430 BOX *
431 ptaGetBoundingRegion(PTA *pta)
432 {
433 l_int32 n, i, x, y, minx, maxx, miny, maxy;
434
435 if (!pta)
436 return (BOX *)ERROR_PTR("pta not defined", __func__, NULL);
437
438 minx = 10000000;
439 miny = 10000000;
440 maxx = -10000000;
441 maxy = -10000000;
442 n = ptaGetCount(pta);
443 for (i = 0; i < n; i++) {
444 ptaGetIPt(pta, i, &x, &y);
445 if (x < minx) minx = x;
446 if (x > maxx) maxx = x;
447 if (y < miny) miny = y;
448 if (y > maxy) maxy = y;
449 }
450
451 return boxCreate(minx, miny, maxx - minx + 1, maxy - miny + 1);
452 }
453
454
455 /*!
456 * \brief ptaGetRange()
457 *
458 * \param[in] pta
459 * \param[out] pminx [optional] min value of x
460 * \param[out] pmaxx [optional] max value of x
461 * \param[out] pminy [optional] min value of y
462 * \param[out] pmaxy [optional] max value of y
463 * \return 0 if OK, 1 on error
464 *
465 * <pre>
466 * Notes:
467 * (1) We can use pts to represent pairs of floating values, that
468 * are not necessarily tied to a two-dimension region. For
469 * example, the pts can represent a general function y(x).
470 * </pre>
471 */
472 l_ok
473 ptaGetRange(PTA *pta,
474 l_float32 *pminx,
475 l_float32 *pmaxx,
476 l_float32 *pminy,
477 l_float32 *pmaxy)
478 {
479 l_int32 n, i;
480 l_float32 x, y, minx, maxx, miny, maxy;
481
482 if (!pminx && !pmaxx && !pminy && !pmaxy)
483 return ERROR_INT("no output requested", __func__, 1);
484 if (pminx) *pminx = 0;
485 if (pmaxx) *pmaxx = 0;
486 if (pminy) *pminy = 0;
487 if (pmaxy) *pmaxy = 0;
488 if (!pta)
489 return ERROR_INT("pta not defined", __func__, 1);
490 if ((n = ptaGetCount(pta)) == 0)
491 return ERROR_INT("no points in pta", __func__, 1);
492
493 ptaGetPt(pta, 0, &x, &y);
494 minx = x;
495 maxx = x;
496 miny = y;
497 maxy = y;
498 for (i = 1; i < n; i++) {
499 ptaGetPt(pta, i, &x, &y);
500 if (x < minx) minx = x;
501 if (x > maxx) maxx = x;
502 if (y < miny) miny = y;
503 if (y > maxy) maxy = y;
504 }
505 if (pminx) *pminx = minx;
506 if (pmaxx) *pmaxx = maxx;
507 if (pminy) *pminy = miny;
508 if (pmaxy) *pmaxy = maxy;
509 return 0;
510 }
511
512
513 /*!
514 * \brief ptaGetInsideBox()
515 *
516 * \param[in] ptas input pts
517 * \param[in] box
518 * \return ptad of pts in ptas that are inside the box, or NULL on error
519 */
520 PTA *
521 ptaGetInsideBox(PTA *ptas,
522 BOX *box)
523 {
524 PTA *ptad;
525 l_int32 n, i, contains;
526 l_float32 x, y;
527
528 if (!ptas)
529 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
530 if (!box)
531 return (PTA *)ERROR_PTR("box not defined", __func__, NULL);
532
533 n = ptaGetCount(ptas);
534 ptad = ptaCreate(0);
535 for (i = 0; i < n; i++) {
536 ptaGetPt(ptas, i, &x, &y);
537 boxContainsPt(box, x, y, &contains);
538 if (contains)
539 ptaAddPt(ptad, x, y);
540 }
541
542 return ptad;
543 }
544
545
546 /*!
547 * \brief pixFindCornerPixels()
548 *
549 * \param[in] pixs 1 bpp
550 * \return pta, or NULL on error
551 *
552 * <pre>
553 * Notes:
554 * (1) Finds the 4 corner-most pixels, as defined by a search
555 * inward from each corner, using a 45 degree line.
556 * </pre>
557 */
558 PTA *
559 pixFindCornerPixels(PIX *pixs)
560 {
561 l_int32 i, j, x, y, w, h, wpl, mindim, found;
562 l_uint32 *data, *line;
563 PTA *pta;
564
565 if (!pixs)
566 return (PTA *)ERROR_PTR("pixs not defined", __func__, NULL);
567 if (pixGetDepth(pixs) != 1)
568 return (PTA *)ERROR_PTR("pixs not 1 bpp", __func__, NULL);
569
570 w = pixGetWidth(pixs);
571 h = pixGetHeight(pixs);
572 mindim = L_MIN(w, h);
573 data = pixGetData(pixs);
574 wpl = pixGetWpl(pixs);
575
576 if ((pta = ptaCreate(4)) == NULL)
577 return (PTA *)ERROR_PTR("pta not made", __func__, NULL);
578
579 for (found = FALSE, i = 0; i < mindim; i++) {
580 for (j = 0; j <= i; j++) {
581 y = i - j;
582 line = data + y * wpl;
583 if (GET_DATA_BIT(line, j)) {
584 ptaAddPt(pta, j, y);
585 found = TRUE;
586 break;
587 }
588 }
589 if (found == TRUE)
590 break;
591 }
592
593 for (found = FALSE, i = 0; i < mindim; i++) {
594 for (j = 0; j <= i; j++) {
595 y = i - j;
596 line = data + y * wpl;
597 x = w - 1 - j;
598 if (GET_DATA_BIT(line, x)) {
599 ptaAddPt(pta, x, y);
600 found = TRUE;
601 break;
602 }
603 }
604 if (found == TRUE)
605 break;
606 }
607
608 for (found = FALSE, i = 0; i < mindim; i++) {
609 for (j = 0; j <= i; j++) {
610 y = h - 1 - i + j;
611 line = data + y * wpl;
612 if (GET_DATA_BIT(line, j)) {
613 ptaAddPt(pta, j, y);
614 found = TRUE;
615 break;
616 }
617 }
618 if (found == TRUE)
619 break;
620 }
621
622 for (found = FALSE, i = 0; i < mindim; i++) {
623 for (j = 0; j <= i; j++) {
624 y = h - 1 - i + j;
625 line = data + y * wpl;
626 x = w - 1 - j;
627 if (GET_DATA_BIT(line, x)) {
628 ptaAddPt(pta, x, y);
629 found = TRUE;
630 break;
631 }
632 }
633 if (found == TRUE)
634 break;
635 }
636
637 return pta;
638 }
639
640
641 /*!
642 * \brief ptaContainsPt()
643 *
644 * \param[in] pta
645 * \param[in] x, y point
646 * \return 1 if contained, 0 otherwise or on error
647 */
648 l_int32
649 ptaContainsPt(PTA *pta,
650 l_int32 x,
651 l_int32 y)
652 {
653 l_int32 i, n, ix, iy;
654
655 if (!pta)
656 return ERROR_INT("pta not defined", __func__, 0);
657
658 n = ptaGetCount(pta);
659 for (i = 0; i < n; i++) {
660 ptaGetIPt(pta, i, &ix, &iy);
661 if (x == ix && y == iy)
662 return 1;
663 }
664 return 0;
665 }
666
667
668 /*!
669 * \brief ptaTestIntersection()
670 *
671 * \param[in] pta1, pta2
672 * \return bval which is 1 if they have any elements in common;
673 * 0 otherwise or on error.
674 */
675 l_int32
676 ptaTestIntersection(PTA *pta1,
677 PTA *pta2)
678 {
679 l_int32 i, j, n1, n2, x1, y1, x2, y2;
680
681 if (!pta1)
682 return ERROR_INT("pta1 not defined", __func__, 0);
683 if (!pta2)
684 return ERROR_INT("pta2 not defined", __func__, 0);
685
686 n1 = ptaGetCount(pta1);
687 n2 = ptaGetCount(pta2);
688 for (i = 0; i < n1; i++) {
689 ptaGetIPt(pta1, i, &x1, &y1);
690 for (j = 0; j < n2; j++) {
691 ptaGetIPt(pta2, i, &x2, &y2);
692 if (x1 == x2 && y1 == y2)
693 return 1;
694 }
695 }
696
697 return 0;
698 }
699
700
701 /*!
702 * \brief ptaTransform()
703 *
704 * \param[in] ptas
705 * \param[in] shiftx, shifty
706 * \param[in] scalex, scaley
707 * \return pta, or NULL on error
708 *
709 * <pre>
710 * Notes:
711 * (1) Shift first, then scale.
712 * </pre>
713 */
714 PTA *
715 ptaTransform(PTA *ptas,
716 l_int32 shiftx,
717 l_int32 shifty,
718 l_float32 scalex,
719 l_float32 scaley)
720 {
721 l_int32 n, i, x, y;
722 PTA *ptad;
723
724 if (!ptas)
725 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
726 n = ptaGetCount(ptas);
727 ptad = ptaCreate(n);
728 for (i = 0; i < n; i++) {
729 ptaGetIPt(ptas, i, &x, &y);
730 x = (l_int32)(scalex * (x + shiftx) + 0.5);
731 y = (l_int32)(scaley * (y + shifty) + 0.5);
732 ptaAddPt(ptad, x, y);
733 }
734
735 return ptad;
736 }
737
738
739 /*!
740 * \brief ptaPtInsidePolygon()
741 *
742 * \param[in] pta vertices of a polygon
743 * \param[in] x, y point to be tested
744 * \param[out] pinside 1 if inside; 0 if outside or on boundary
745 * \return 1 if OK, 0 on error
746 *
747 * The abs value of the sum of the angles subtended from a point by
748 * the sides of a polygon, when taken in order traversing the polygon,
749 * is 0 if the point is outside the polygon and 2*pi if inside.
750 * The sign will be positive if traversed cw and negative if ccw.
751 */
752 l_int32
753 ptaPtInsidePolygon(PTA *pta,
754 l_float32 x,
755 l_float32 y,
756 l_int32 *pinside)
757 {
758 l_int32 i, n;
759 l_float32 sum, x1, y1, x2, y2, xp1, yp1, xp2, yp2;
760
761 if (!pinside)
762 return ERROR_INT("&inside not defined", __func__, 1);
763 *pinside = 0;
764 if (!pta)
765 return ERROR_INT("pta not defined", __func__, 1);
766
767 /* Think of (x1,y1) as the end point of a vector that starts
768 * from the origin (0,0), and ditto for (x2,y2). */
769 n = ptaGetCount(pta);
770 sum = 0.0;
771 for (i = 0; i < n; i++) {
772 ptaGetPt(pta, i, &xp1, &yp1);
773 ptaGetPt(pta, (i + 1) % n, &xp2, &yp2);
774 x1 = xp1 - x;
775 y1 = yp1 - y;
776 x2 = xp2 - x;
777 y2 = yp2 - y;
778 sum += l_angleBetweenVectors(x1, y1, x2, y2);
779 }
780
781 if (L_ABS(sum) > M_PI)
782 *pinside = 1;
783 return 0;
784 }
785
786
787 /*!
788 * \brief l_angleBetweenVectors()
789 *
790 * \param[in] x1, y1 end point of first vector
791 * \param[in] x2, y2 end point of second vector
792 * \return angle radians, or 0.0 on error
793 *
794 * <pre>
795 * Notes:
796 * (1) This gives the angle between two vectors, going between
797 * vector1 (x1,y1) and vector2 (x2,y2). The angle is swept
798 * out from 1 --> 2. If this is clockwise, the angle is
799 * positive, but the result is folded into the interval [-pi, pi].
800 * </pre>
801 */
802 l_float32
803 l_angleBetweenVectors(l_float32 x1,
804 l_float32 y1,
805 l_float32 x2,
806 l_float32 y2)
807 {
808 l_float64 ang;
809
810 ang = atan2(y2, x2) - atan2(y1, x1);
811 if (ang > M_PI) ang -= 2.0 * M_PI;
812 if (ang < -M_PI) ang += 2.0 * M_PI;
813 return ang;
814 }
815
816
817 /*!
818 * \brief ptaPolygonIsConvex()
819 *
820 * \param[in] pta corners of polygon
821 * \param[out] pisconvex 1 if convex; 0 otherwise
822 * \return 0 if OK, 1 on error
823 *
824 * <pre>
825 * Notes:
826 * (1) A Pta of size n describes a polygon with n sides, where
827 * the n-th side goes from point[n - 1] to point[0].
828 * (2) The pta must describe a CLOCKWISE traversal of the boundary
829 * of the polygon.
830 * (3) Algorithm: traversing the boundary in a cw direction, the
831 * polygon interior is always on the right. If the polygon is
832 * convex, for each set of 3 points, the third point is either
833 * on the ray extending from the first point and going through
834 * the second point, or to the right of it.
835 * </pre>
836 */
837 l_int32
838 ptaPolygonIsConvex(PTA *pta,
839 l_int32 *pisconvex)
840 {
841 l_int32 i, n;
842 l_float32 x0, y0, x1, y1, x2, y2;
843 l_float64 cprod;
844
845 if (!pisconvex)
846 return ERROR_INT("&isconvex not defined", __func__, 1);
847 *pisconvex = 0;
848 if (!pta)
849 return ERROR_INT("pta not defined", __func__, 1);
850 if ((n = ptaGetCount(pta)) < 3)
851 return ERROR_INT("pta has < 3 pts", __func__, 1);
852
853 for (i = 0; i < n; i++) {
854 ptaGetPt(pta, i, &x0, &y0);
855 ptaGetPt(pta, (i + 1) % n, &x1, &y1);
856 ptaGetPt(pta, (i + 2) % n, &x2, &y2);
857 /* The vector v02 from p0 to p2 must be to the right of the
858 vector v01 from p0 to p1. This is true if the cross
859 product v02 x v01 > 0. In coordinates:
860 v02x * v01y - v01x * v02y, where
861 v01x = x1 - x0, v01y = y1 - y0,
862 v02x = x2 - x0, v02y = y2 - y0 */
863 cprod = (x2 - x0) * (y1 - y0) - (x1 - x0) * (y2 - y0);
864 if (cprod < -0.0001) /* small delta for float accuracy; test fails */
865 return 0;
866 }
867 *pisconvex = 1;
868 return 0;
869 }
870
871
872 /*---------------------------------------------------------------------*
873 * Min/max and filtering *
874 *---------------------------------------------------------------------*/
875 /*!
876 * \brief ptaGetMinMax()
877 *
878 * \param[in] pta
879 * \param[out] pxmin [optional] min of x
880 * \param[out] pymin [optional] min of y
881 * \param[out] pxmax [optional] max of x
882 * \param[out] pymax [optional] max of y
883 * \return 0 if OK, 1 on error. If pta is empty, requested
884 * values are returned as -1.0.
885 */
886 l_ok
887 ptaGetMinMax(PTA *pta,
888 l_float32 *pxmin,
889 l_float32 *pymin,
890 l_float32 *pxmax,
891 l_float32 *pymax)
892 {
893 l_int32 i, n;
894 l_float32 x, y, xmin, ymin, xmax, ymax;
895
896 if (pxmin) *pxmin = -1.0;
897 if (pymin) *pymin = -1.0;
898 if (pxmax) *pxmax = -1.0;
899 if (pymax) *pymax = -1.0;
900 if (!pta)
901 return ERROR_INT("pta not defined", __func__, 1);
902 if (!pxmin && !pxmax && !pymin && !pymax)
903 return ERROR_INT("no output requested", __func__, 1);
904 if ((n = ptaGetCount(pta)) == 0) {
905 L_WARNING("pta is empty\n", __func__);
906 return 0;
907 }
908
909 xmin = ymin = 1.0e20f;
910 xmax = ymax = -1.0e20f;
911 for (i = 0; i < n; i++) {
912 ptaGetPt(pta, i, &x, &y);
913 if (x < xmin) xmin = x;
914 if (y < ymin) ymin = y;
915 if (x > xmax) xmax = x;
916 if (y > ymax) ymax = y;
917 }
918 if (pxmin) *pxmin = xmin;
919 if (pymin) *pymin = ymin;
920 if (pxmax) *pxmax = xmax;
921 if (pymax) *pymax = ymax;
922 return 0;
923 }
924
925
926 /*!
927 * \brief ptaSelectByValue()
928 *
929 * \param[in] ptas
930 * \param[in] xth, yth threshold values
931 * \param[in] type L_SELECT_XVAL, L_SELECT_YVAL,
932 * L_SELECT_IF_EITHER, L_SELECT_IF_BOTH
933 * \param[in] relation L_SELECT_IF_LT, L_SELECT_IF_GT,
934 * L_SELECT_IF_LTE, L_SELECT_IF_GTE
935 * \return ptad filtered set, or NULL on error
936 */
937 PTA *
938 ptaSelectByValue(PTA *ptas,
939 l_float32 xth,
940 l_float32 yth,
941 l_int32 type,
942 l_int32 relation)
943 {
944 l_int32 i, n;
945 l_float32 x, y;
946 PTA *ptad;
947
948 if (!ptas)
949 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
950 if (ptaGetCount(ptas) == 0) {
951 L_WARNING("ptas is empty\n", __func__);
952 return ptaCopy(ptas);
953 }
954 if (type != L_SELECT_XVAL && type != L_SELECT_YVAL &&
955 type != L_SELECT_IF_EITHER && type != L_SELECT_IF_BOTH)
956 return (PTA *)ERROR_PTR("invalid type", __func__, NULL);
957 if (relation != L_SELECT_IF_LT && relation != L_SELECT_IF_GT &&
958 relation != L_SELECT_IF_LTE && relation != L_SELECT_IF_GTE)
959 return (PTA *)ERROR_PTR("invalid relation", __func__, NULL);
960
961 n = ptaGetCount(ptas);
962 ptad = ptaCreate(n);
963 for (i = 0; i < n; i++) {
964 ptaGetPt(ptas, i, &x, &y);
965 if (type == L_SELECT_XVAL) {
966 if ((relation == L_SELECT_IF_LT && x < xth) ||
967 (relation == L_SELECT_IF_GT && x > xth) ||
968 (relation == L_SELECT_IF_LTE && x <= xth) ||
969 (relation == L_SELECT_IF_GTE && x >= xth))
970 ptaAddPt(ptad, x, y);
971 } else if (type == L_SELECT_YVAL) {
972 if ((relation == L_SELECT_IF_LT && y < yth) ||
973 (relation == L_SELECT_IF_GT && y > yth) ||
974 (relation == L_SELECT_IF_LTE && y <= yth) ||
975 (relation == L_SELECT_IF_GTE && y >= yth))
976 ptaAddPt(ptad, x, y);
977 } else if (type == L_SELECT_IF_EITHER) {
978 if (((relation == L_SELECT_IF_LT) && (x < xth || y < yth)) ||
979 ((relation == L_SELECT_IF_GT) && (x > xth || y > yth)) ||
980 ((relation == L_SELECT_IF_LTE) && (x <= xth || y <= yth)) ||
981 ((relation == L_SELECT_IF_GTE) && (x >= xth || y >= yth)))
982 ptaAddPt(ptad, x, y);
983 } else { /* L_SELECT_IF_BOTH */
984 if (((relation == L_SELECT_IF_LT) && (x < xth && y < yth)) ||
985 ((relation == L_SELECT_IF_GT) && (x > xth && y > yth)) ||
986 ((relation == L_SELECT_IF_LTE) && (x <= xth && y <= yth)) ||
987 ((relation == L_SELECT_IF_GTE) && (x >= xth && y >= yth)))
988 ptaAddPt(ptad, x, y);
989 }
990 }
991
992 return ptad;
993 }
994
995
996 /*!
997 * \brief ptaCropToMask()
998 *
999 * \param[in] ptas input pta
1000 * \param[in] pixm 1 bpp mask
1001 * \return ptad with only pts under the mask fg, or NULL on error
1002 */
1003 PTA *
1004 ptaCropToMask(PTA *ptas,
1005 PIX *pixm)
1006 {
1007 l_int32 i, n, x, y;
1008 l_uint32 val;
1009 PTA *ptad;
1010
1011 if (!ptas)
1012 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
1013 if (!pixm || pixGetDepth(pixm) != 1)
1014 return (PTA *)ERROR_PTR("pixm undefined or not 1 bpp", __func__, NULL);
1015 if (ptaGetCount(ptas) == 0) {
1016 L_INFO("ptas is empty\n", __func__);
1017 return ptaCopy(ptas);
1018 }
1019
1020 n = ptaGetCount(ptas);
1021 ptad = ptaCreate(n);
1022 for (i = 0; i < n; i++) {
1023 ptaGetIPt(ptas, i, &x, &y);
1024 pixGetPixel(pixm, x, y, &val);
1025 if (val == 1)
1026 ptaAddPt(ptad, x, y);
1027 }
1028 return ptad;
1029 }
1030
1031
1032 /*---------------------------------------------------------------------*
1033 * Least Squares Fit *
1034 *---------------------------------------------------------------------*/
1035 /*!
1036 * \brief ptaGetLinearLSF()
1037 *
1038 * \param[in] pta
1039 * \param[out] pa [optional] slope a of least square fit: y = ax + b
1040 * \param[out] pb [optional] intercept b of least square fit
1041 * \param[out] pnafit [optional] numa of least square fit
1042 * \return 0 if OK, 1 on error
1043 *
1044 * <pre>
1045 * Notes:
1046 * (1) Either or both &a and &b must be input. They determine the
1047 * type of line that is fit.
1048 * (2) If both &a and &b are defined, this returns a and b that minimize:
1049 *
1050 * sum (yi - axi -b)^2
1051 * i
1052 *
1053 * The method is simple: differentiate this expression w/rt a and b,
1054 * and solve the resulting two equations for a and b in terms of
1055 * various sums over the input data (xi, yi).
1056 * (3) We also allow two special cases, where either a = 0 or b = 0:
1057 * (a) If &a is given and &b = null, find the linear LSF that
1058 * goes through the origin (b = 0).
1059 * (b) If &b is given and &a = null, find the linear LSF with
1060 * zero slope (a = 0).
1061 * (4) If &nafit is defined, this returns an array of fitted values,
1062 * corresponding to the two implicit Numa arrays (nax and nay) in pta.
1063 * Thus, just as you can plot the data in pta as nay vs. nax,
1064 * you can plot the linear least square fit as nafit vs. nax.
1065 * Get the nax array using ptaGetArrays(pta, &nax, NULL);
1066 * </pre>
1067 */
1068 l_ok
1069 ptaGetLinearLSF(PTA *pta,
1070 l_float32 *pa,
1071 l_float32 *pb,
1072 NUMA **pnafit)
1073 {
1074 l_int32 n, i;
1075 l_float32 a, b, factor, sx, sy, sxx, sxy, val;
1076 l_float32 *xa, *ya;
1077
1078 if (pa) *pa = 0.0;
1079 if (pb) *pb = 0.0;
1080 if (pnafit) *pnafit = NULL;
1081 if (!pa && !pb && !pnafit)
1082 return ERROR_INT("no output requested", __func__, 1);
1083 if (!pta)
1084 return ERROR_INT("pta not defined", __func__, 1);
1085 if ((n = ptaGetCount(pta)) < 2)
1086 return ERROR_INT("less than 2 pts found", __func__, 1);
1087
1088 xa = pta->x; /* not a copy */
1089 ya = pta->y; /* not a copy */
1090 sx = sy = sxx = sxy = 0.;
1091 if (pa && pb) { /* general line */
1092 for (i = 0; i < n; i++) {
1093 sx += xa[i];
1094 sy += ya[i];
1095 sxx += xa[i] * xa[i];
1096 sxy += xa[i] * ya[i];
1097 }
1098 factor = n * sxx - sx * sx;
1099 if (factor == 0.0)
1100 return ERROR_INT("no solution found", __func__, 1);
1101 factor = 1.f / factor;
1102
1103 a = factor * ((l_float32)n * sxy - sx * sy);
1104 b = factor * (sxx * sy - sx * sxy);
1105 } else if (pa) { /* b = 0; line through origin */
1106 for (i = 0; i < n; i++) {
1107 sxx += xa[i] * xa[i];
1108 sxy += xa[i] * ya[i];
1109 }
1110 if (sxx == 0.0)
1111 return ERROR_INT("no solution found", __func__, 1);
1112 a = sxy / sxx;
1113 b = 0.0;
1114 } else { /* a = 0; horizontal line */
1115 for (i = 0; i < n; i++)
1116 sy += ya[i];
1117 a = 0.0;
1118 b = sy / (l_float32)n;
1119 }
1120
1121 if (pnafit) {
1122 *pnafit = numaCreate(n);
1123 for (i = 0; i < n; i++) {
1124 val = a * xa[i] + b;
1125 numaAddNumber(*pnafit, val);
1126 }
1127 }
1128
1129 if (pa) *pa = a;
1130 if (pb) *pb = b;
1131 return 0;
1132 }
1133
1134
1135 /*!
1136 * \brief ptaGetQuadraticLSF()
1137 *
1138 * \param[in] pta
1139 * \param[out] pa [optional] coeff a of LSF: y = ax^2 + bx + c
1140 * \param[out] pb [optional] coeff b of LSF: y = ax^2 + bx + c
1141 * \param[out] pc [optional] coeff c of LSF: y = ax^2 + bx + c
1142 * \param[out] pnafit [optional] numa of least square fit
1143 * \return 0 if OK, 1 on error
1144 *
1145 * <pre>
1146 * Notes:
1147 * (1) This does a quadratic least square fit to the set of points
1148 * in %pta. That is, it finds coefficients a, b and c that minimize:
1149 *
1150 * sum (yi - a*xi*xi -b*xi -c)^2
1151 * i
1152 *
1153 * The method is simple: differentiate this expression w/rt
1154 * a, b and c, and solve the resulting three equations for these
1155 * coefficients in terms of various sums over the input data (xi, yi).
1156 * The three equations are in the form:
1157 * f[0][0]a + f[0][1]b + f[0][2]c = g[0]
1158 * f[1][0]a + f[1][1]b + f[1][2]c = g[1]
1159 * f[2][0]a + f[2][1]b + f[2][2]c = g[2]
1160 * (2) If &nafit is defined, this returns an array of fitted values,
1161 * corresponding to the two implicit Numa arrays (nax and nay) in pta.
1162 * Thus, just as you can plot the data in pta as nay vs. nax,
1163 * you can plot the linear least square fit as nafit vs. nax.
1164 * Get the nax array using ptaGetArrays(pta, &nax, NULL);
1165 * </pre>
1166 */
1167 l_ok
1168 ptaGetQuadraticLSF(PTA *pta,
1169 l_float32 *pa,
1170 l_float32 *pb,
1171 l_float32 *pc,
1172 NUMA **pnafit)
1173 {
1174 l_int32 n, i, ret;
1175 l_float32 x, y, sx, sy, sx2, sx3, sx4, sxy, sx2y;
1176 l_float32 *xa, *ya;
1177 l_float32 *f[3];
1178 l_float32 g[3];
1179
1180 if (pa) *pa = 0.0;
1181 if (pb) *pb = 0.0;
1182 if (pc) *pc = 0.0;
1183 if (pnafit) *pnafit = NULL;
1184 if (!pa && !pb && !pc && !pnafit)
1185 return ERROR_INT("no output requested", __func__, 1);
1186 if (!pta)
1187 return ERROR_INT("pta not defined", __func__, 1);
1188 if ((n = ptaGetCount(pta)) < 3)
1189 return ERROR_INT("less than 3 pts found", __func__, 1);
1190
1191 xa = pta->x; /* not a copy */
1192 ya = pta->y; /* not a copy */
1193 sx = sy = sx2 = sx3 = sx4 = sxy = sx2y = 0.;
1194 for (i = 0; i < n; i++) {
1195 x = xa[i];
1196 y = ya[i];
1197 sx += x;
1198 sy += y;
1199 sx2 += x * x;
1200 sx3 += x * x * x;
1201 sx4 += x * x * x * x;
1202 sxy += x * y;
1203 sx2y += x * x * y;
1204 }
1205
1206 for (i = 0; i < 3; i++)
1207 f[i] = (l_float32 *)LEPT_CALLOC(3, sizeof(l_float32));
1208 f[0][0] = sx4;
1209 f[0][1] = sx3;
1210 f[0][2] = sx2;
1211 f[1][0] = sx3;
1212 f[1][1] = sx2;
1213 f[1][2] = sx;
1214 f[2][0] = sx2;
1215 f[2][1] = sx;
1216 f[2][2] = n;
1217 g[0] = sx2y;
1218 g[1] = sxy;
1219 g[2] = sy;
1220
1221 /* Solve for the unknowns, also putting f-inverse into f */
1222 ret = gaussjordan(f, g, 3);
1223 for (i = 0; i < 3; i++)
1224 LEPT_FREE(f[i]);
1225 if (ret)
1226 return ERROR_INT("quadratic solution failed", __func__, 1);
1227
1228 if (pa) *pa = g[0];
1229 if (pb) *pb = g[1];
1230 if (pc) *pc = g[2];
1231 if (pnafit) {
1232 *pnafit = numaCreate(n);
1233 for (i = 0; i < n; i++) {
1234 x = xa[i];
1235 y = g[0] * x * x + g[1] * x + g[2];
1236 numaAddNumber(*pnafit, y);
1237 }
1238 }
1239 return 0;
1240 }
1241
1242
1243 /*!
1244 * \brief ptaGetCubicLSF()
1245 *
1246 * \param[in] pta
1247 * \param[out] pa [optional] coeff a of LSF: y = ax^3 + bx^2 + cx + d
1248 * \param[out] pb [optional] coeff b of LSF
1249 * \param[out] pc [optional] coeff c of LSF
1250 * \param[out] pd [optional] coeff d of LSF
1251 * \param[out] pnafit [optional] numa of least square fit
1252 * \return 0 if OK, 1 on error
1253 *
1254 * <pre>
1255 * Notes:
1256 * (1) This does a cubic least square fit to the set of points
1257 * in %pta. That is, it finds coefficients a, b, c and d
1258 * that minimize:
1259 *
1260 * sum (yi - a*xi*xi*xi -b*xi*xi -c*xi - d)^2
1261 * i
1262 *
1263 * Differentiate this expression w/rt a, b, c and d, and solve
1264 * the resulting four equations for these coefficients in
1265 * terms of various sums over the input data (xi, yi).
1266 * The four equations are in the form:
1267 * f[0][0]a + f[0][1]b + f[0][2]c + f[0][3] = g[0]
1268 * f[1][0]a + f[1][1]b + f[1][2]c + f[1][3] = g[1]
1269 * f[2][0]a + f[2][1]b + f[2][2]c + f[2][3] = g[2]
1270 * f[3][0]a + f[3][1]b + f[3][2]c + f[3][3] = g[3]
1271 * (2) If &nafit is defined, this returns an array of fitted values,
1272 * corresponding to the two implicit Numa arrays (nax and nay) in pta.
1273 * Thus, just as you can plot the data in pta as nay vs. nax,
1274 * you can plot the linear least square fit as nafit vs. nax.
1275 * Get the nax array using ptaGetArrays(pta, &nax, NULL);
1276 * </pre>
1277 */
1278 l_ok
1279 ptaGetCubicLSF(PTA *pta,
1280 l_float32 *pa,
1281 l_float32 *pb,
1282 l_float32 *pc,
1283 l_float32 *pd,
1284 NUMA **pnafit)
1285 {
1286 l_int32 n, i, ret;
1287 l_float32 x, y, sx, sy, sx2, sx3, sx4, sx5, sx6, sxy, sx2y, sx3y;
1288 l_float32 *xa, *ya;
1289 l_float32 *f[4];
1290 l_float32 g[4];
1291
1292 if (pa) *pa = 0.0;
1293 if (pb) *pb = 0.0;
1294 if (pc) *pc = 0.0;
1295 if (pd) *pd = 0.0;
1296 if (pnafit) *pnafit = NULL;
1297 if (!pa && !pb && !pc && !pd && !pnafit)
1298 return ERROR_INT("no output requested", __func__, 1);
1299 if (!pta)
1300 return ERROR_INT("pta not defined", __func__, 1);
1301 if ((n = ptaGetCount(pta)) < 4)
1302 return ERROR_INT("less than 4 pts found", __func__, 1);
1303
1304 xa = pta->x; /* not a copy */
1305 ya = pta->y; /* not a copy */
1306 sx = sy = sx2 = sx3 = sx4 = sx5 = sx6 = sxy = sx2y = sx3y = 0.;
1307 for (i = 0; i < n; i++) {
1308 x = xa[i];
1309 y = ya[i];
1310 sx += x;
1311 sy += y;
1312 sx2 += x * x;
1313 sx3 += x * x * x;
1314 sx4 += x * x * x * x;
1315 sx5 += x * x * x * x * x;
1316 sx6 += x * x * x * x * x * x;
1317 sxy += x * y;
1318 sx2y += x * x * y;
1319 sx3y += x * x * x * y;
1320 }
1321
1322 for (i = 0; i < 4; i++)
1323 f[i] = (l_float32 *)LEPT_CALLOC(4, sizeof(l_float32));
1324 f[0][0] = sx6;
1325 f[0][1] = sx5;
1326 f[0][2] = sx4;
1327 f[0][3] = sx3;
1328 f[1][0] = sx5;
1329 f[1][1] = sx4;
1330 f[1][2] = sx3;
1331 f[1][3] = sx2;
1332 f[2][0] = sx4;
1333 f[2][1] = sx3;
1334 f[2][2] = sx2;
1335 f[2][3] = sx;
1336 f[3][0] = sx3;
1337 f[3][1] = sx2;
1338 f[3][2] = sx;
1339 f[3][3] = n;
1340 g[0] = sx3y;
1341 g[1] = sx2y;
1342 g[2] = sxy;
1343 g[3] = sy;
1344
1345 /* Solve for the unknowns, also putting f-inverse into f */
1346 ret = gaussjordan(f, g, 4);
1347 for (i = 0; i < 4; i++)
1348 LEPT_FREE(f[i]);
1349 if (ret)
1350 return ERROR_INT("cubic solution failed", __func__, 1);
1351
1352 if (pa) *pa = g[0];
1353 if (pb) *pb = g[1];
1354 if (pc) *pc = g[2];
1355 if (pd) *pd = g[3];
1356 if (pnafit) {
1357 *pnafit = numaCreate(n);
1358 for (i = 0; i < n; i++) {
1359 x = xa[i];
1360 y = g[0] * x * x * x + g[1] * x * x + g[2] * x + g[3];
1361 numaAddNumber(*pnafit, y);
1362 }
1363 }
1364 return 0;
1365 }
1366
1367
1368 /*!
1369 * \brief ptaGetQuarticLSF()
1370 *
1371 * \param[in] pta
1372 * \param[out] pa [optional] coeff a of LSF:
1373 * y = ax^4 + bx^3 + cx^2 + dx + e
1374 * \param[out] pb [optional] coeff b of LSF
1375 * \param[out] pc [optional] coeff c of LSF
1376 * \param[out] pd [optional] coeff d of LSF
1377 * \param[out] pe [optional] coeff e of LSF
1378 * \param[out] pnafit [optional] numa of least square fit
1379 * \return 0 if OK, 1 on error
1380 *
1381 * <pre>
1382 * Notes:
1383 * (1) This does a quartic least square fit to the set of points
1384 * in %pta. That is, it finds coefficients a, b, c, d and 3
1385 * that minimize:
1386 *
1387 * sum (yi - a*xi*xi*xi*xi -b*xi*xi*xi -c*xi*xi - d*xi - e)^2
1388 * i
1389 *
1390 * Differentiate this expression w/rt a, b, c, d and e, and solve
1391 * the resulting five equations for these coefficients in
1392 * terms of various sums over the input data (xi, yi).
1393 * The five equations are in the form:
1394 * f[0][0]a + f[0][1]b + f[0][2]c + f[0][3] + f[0][4] = g[0]
1395 * f[1][0]a + f[1][1]b + f[1][2]c + f[1][3] + f[1][4] = g[1]
1396 * f[2][0]a + f[2][1]b + f[2][2]c + f[2][3] + f[2][4] = g[2]
1397 * f[3][0]a + f[3][1]b + f[3][2]c + f[3][3] + f[3][4] = g[3]
1398 * f[4][0]a + f[4][1]b + f[4][2]c + f[4][3] + f[4][4] = g[4]
1399 * (2) If &nafit is defined, this returns an array of fitted values,
1400 * corresponding to the two implicit Numa arrays (nax and nay) in pta.
1401 * Thus, just as you can plot the data in pta as nay vs. nax,
1402 * you can plot the linear least square fit as nafit vs. nax.
1403 * Get the nax array using ptaGetArrays(pta, &nax, NULL);
1404 * </pre>
1405 */
1406 l_ok
1407 ptaGetQuarticLSF(PTA *pta,
1408 l_float32 *pa,
1409 l_float32 *pb,
1410 l_float32 *pc,
1411 l_float32 *pd,
1412 l_float32 *pe,
1413 NUMA **pnafit)
1414 {
1415 l_int32 n, i, ret;
1416 l_float32 x, y, sx, sy, sx2, sx3, sx4, sx5, sx6, sx7, sx8;
1417 l_float32 sxy, sx2y, sx3y, sx4y;
1418 l_float32 *xa, *ya;
1419 l_float32 *f[5];
1420 l_float32 g[5];
1421
1422 if (pa) *pa = 0.0;
1423 if (pb) *pb = 0.0;
1424 if (pc) *pc = 0.0;
1425 if (pd) *pd = 0.0;
1426 if (pe) *pe = 0.0;
1427 if (pnafit) *pnafit = NULL;
1428 if (!pa && !pb && !pc && !pd && !pe && !pnafit)
1429 return ERROR_INT("no output requested", __func__, 1);
1430 if (!pta)
1431 return ERROR_INT("pta not defined", __func__, 1);
1432 if ((n = ptaGetCount(pta)) < 5)
1433 return ERROR_INT("less than 5 pts found", __func__, 1);
1434
1435 xa = pta->x; /* not a copy */
1436 ya = pta->y; /* not a copy */
1437 sx = sy = sx2 = sx3 = sx4 = sx5 = sx6 = sx7 = sx8 = 0;
1438 sxy = sx2y = sx3y = sx4y = 0.;
1439 for (i = 0; i < n; i++) {
1440 x = xa[i];
1441 y = ya[i];
1442 sx += x;
1443 sy += y;
1444 sx2 += x * x;
1445 sx3 += x * x * x;
1446 sx4 += x * x * x * x;
1447 sx5 += x * x * x * x * x;
1448 sx6 += x * x * x * x * x * x;
1449 sx7 += x * x * x * x * x * x * x;
1450 sx8 += x * x * x * x * x * x * x * x;
1451 sxy += x * y;
1452 sx2y += x * x * y;
1453 sx3y += x * x * x * y;
1454 sx4y += x * x * x * x * y;
1455 }
1456
1457 for (i = 0; i < 5; i++)
1458 f[i] = (l_float32 *)LEPT_CALLOC(5, sizeof(l_float32));
1459 f[0][0] = sx8;
1460 f[0][1] = sx7;
1461 f[0][2] = sx6;
1462 f[0][3] = sx5;
1463 f[0][4] = sx4;
1464 f[1][0] = sx7;
1465 f[1][1] = sx6;
1466 f[1][2] = sx5;
1467 f[1][3] = sx4;
1468 f[1][4] = sx3;
1469 f[2][0] = sx6;
1470 f[2][1] = sx5;
1471 f[2][2] = sx4;
1472 f[2][3] = sx3;
1473 f[2][4] = sx2;
1474 f[3][0] = sx5;
1475 f[3][1] = sx4;
1476 f[3][2] = sx3;
1477 f[3][3] = sx2;
1478 f[3][4] = sx;
1479 f[4][0] = sx4;
1480 f[4][1] = sx3;
1481 f[4][2] = sx2;
1482 f[4][3] = sx;
1483 f[4][4] = n;
1484 g[0] = sx4y;
1485 g[1] = sx3y;
1486 g[2] = sx2y;
1487 g[3] = sxy;
1488 g[4] = sy;
1489
1490 /* Solve for the unknowns, also putting f-inverse into f */
1491 ret = gaussjordan(f, g, 5);
1492 for (i = 0; i < 5; i++)
1493 LEPT_FREE(f[i]);
1494 if (ret)
1495 return ERROR_INT("quartic solution failed", __func__, 1);
1496
1497 if (pa) *pa = g[0];
1498 if (pb) *pb = g[1];
1499 if (pc) *pc = g[2];
1500 if (pd) *pd = g[3];
1501 if (pe) *pe = g[4];
1502 if (pnafit) {
1503 *pnafit = numaCreate(n);
1504 for (i = 0; i < n; i++) {
1505 x = xa[i];
1506 y = g[0] * x * x * x * x + g[1] * x * x * x + g[2] * x * x
1507 + g[3] * x + g[4];
1508 numaAddNumber(*pnafit, y);
1509 }
1510 }
1511 return 0;
1512 }
1513
1514
1515 /*!
1516 * \brief ptaNoisyLinearLSF()
1517 *
1518 * \param[in] pta
1519 * \param[in] factor reject outliers with error greater than this
1520 * number of medians; typically ~ 3
1521 * \param[out] pptad [optional] with outliers removed
1522 * \param[out] pa [optional] slope a of least square fit: y = ax + b
1523 * \param[out] pb [optional] intercept b of least square fit
1524 * \param[out] pmederr [optional] median error
1525 * \param[out] pnafit [optional] numa of least square fit to ptad
1526 * \return 0 if OK, 1 on error
1527 *
1528 * <pre>
1529 * Notes:
1530 * (1) This does a linear least square fit to the set of points
1531 * in %pta. It then evaluates the errors and removes points
1532 * whose error is >= factor * median_error. It then re-runs
1533 * the linear LSF on the resulting points.
1534 * (2) Either or both &a and &b must be input. They determine the
1535 * type of line that is fit.
1536 * (3) The median error can give an indication of how good the fit
1537 * is likely to be.
1538 * </pre>
1539 */
1540 l_ok
1541 ptaNoisyLinearLSF(PTA *pta,
1542 l_float32 factor,
1543 PTA **pptad,
1544 l_float32 *pa,
1545 l_float32 *pb,
1546 l_float32 *pmederr,
1547 NUMA **pnafit)
1548 {
1549 l_int32 n, i, ret;
1550 l_float32 x, y, yf, val, mederr;
1551 NUMA *nafit, *naerror;
1552 PTA *ptad;
1553
1554 if (pptad) *pptad = NULL;
1555 if (pa) *pa = 0.0;
1556 if (pb) *pb = 0.0;
1557 if (pmederr) *pmederr = 0.0;
1558 if (pnafit) *pnafit = NULL;
1559 if (!pptad && !pa && !pb && !pnafit)
1560 return ERROR_INT("no output requested", __func__, 1);
1561 if (!pta)
1562 return ERROR_INT("pta not defined", __func__, 1);
1563 if (factor <= 0.0)
1564 return ERROR_INT("factor must be > 0.0", __func__, 1);
1565 if ((n = ptaGetCount(pta)) < 3)
1566 return ERROR_INT("less than 2 pts found", __func__, 1);
1567
1568 if (ptaGetLinearLSF(pta, pa, pb, &nafit) != 0)
1569 return ERROR_INT("error in linear LSF", __func__, 1);
1570
1571 /* Get the median error */
1572 naerror = numaCreate(n);
1573 for (i = 0; i < n; i++) {
1574 ptaGetPt(pta, i, &x, &y);
1575 numaGetFValue(nafit, i, &yf);
1576 numaAddNumber(naerror, L_ABS(y - yf));
1577 }
1578 numaGetMedian(naerror, &mederr);
1579 if (pmederr) *pmederr = mederr;
1580 numaDestroy(&nafit);
1581
1582 /* Remove outliers */
1583 ptad = ptaCreate(n);
1584 for (i = 0; i < n; i++) {
1585 ptaGetPt(pta, i, &x, &y);
1586 numaGetFValue(naerror, i, &val);
1587 if (val <= factor * mederr) /* <= in case mederr = 0 */
1588 ptaAddPt(ptad, x, y);
1589 }
1590 numaDestroy(&naerror);
1591
1592 /* Do LSF again */
1593 ret = ptaGetLinearLSF(ptad, pa, pb, pnafit);
1594 if (pptad)
1595 *pptad = ptad;
1596 else
1597 ptaDestroy(&ptad);
1598
1599 return ret;
1600 }
1601
1602
1603 /*!
1604 * \brief ptaNoisyQuadraticLSF()
1605 *
1606 * \param[in] pta
1607 * \param[in] factor reject outliers with error greater than this
1608 * number of medians; typically ~ 3
1609 * \param[out] pptad [optional] with outliers removed
1610 * \param[out] pa [optional] coeff a of LSF: y = ax^2 + bx + c
1611 * \param[out] pb [optional] coeff b of LSF: y = ax^2 + bx + c
1612 * \param[out] pc [optional] coeff c of LSF: y = ax^2 + bx + c
1613 * \param[out] pmederr [optional] median error
1614 * \param[out] pnafit [optional] numa of least square fit to ptad
1615 * \return 0 if OK, 1 on error
1616 *
1617 * <pre>
1618 * Notes:
1619 * (1) This does a quadratic least square fit to the set of points
1620 * in %pta. It then evaluates the errors and removes points
1621 * whose error is >= factor * median_error. It then re-runs
1622 * a quadratic LSF on the resulting points.
1623 * </pre>
1624 */
1625 l_ok
1626 ptaNoisyQuadraticLSF(PTA *pta,
1627 l_float32 factor,
1628 PTA **pptad,
1629 l_float32 *pa,
1630 l_float32 *pb,
1631 l_float32 *pc,
1632 l_float32 *pmederr,
1633 NUMA **pnafit)
1634 {
1635 l_int32 n, i, ret;
1636 l_float32 x, y, yf, val, mederr;
1637 NUMA *nafit, *naerror;
1638 PTA *ptad;
1639
1640 if (pptad) *pptad = NULL;
1641 if (pa) *pa = 0.0;
1642 if (pb) *pb = 0.0;
1643 if (pc) *pc = 0.0;
1644 if (pmederr) *pmederr = 0.0;
1645 if (pnafit) *pnafit = NULL;
1646 if (!pptad && !pa && !pb && !pc && !pnafit)
1647 return ERROR_INT("no output requested", __func__, 1);
1648 if (factor <= 0.0)
1649 return ERROR_INT("factor must be > 0.0", __func__, 1);
1650 if (!pta)
1651 return ERROR_INT("pta not defined", __func__, 1);
1652 if ((n = ptaGetCount(pta)) < 3)
1653 return ERROR_INT("less than 3 pts found", __func__, 1);
1654
1655 if (ptaGetQuadraticLSF(pta, NULL, NULL, NULL, &nafit) != 0)
1656 return ERROR_INT("error in quadratic LSF", __func__, 1);
1657
1658 /* Get the median error */
1659 naerror = numaCreate(n);
1660 for (i = 0; i < n; i++) {
1661 ptaGetPt(pta, i, &x, &y);
1662 numaGetFValue(nafit, i, &yf);
1663 numaAddNumber(naerror, L_ABS(y - yf));
1664 }
1665 numaGetMedian(naerror, &mederr);
1666 if (pmederr) *pmederr = mederr;
1667 numaDestroy(&nafit);
1668
1669 /* Remove outliers */
1670 ptad = ptaCreate(n);
1671 for (i = 0; i < n; i++) {
1672 ptaGetPt(pta, i, &x, &y);
1673 numaGetFValue(naerror, i, &val);
1674 if (val <= factor * mederr) /* <= in case mederr = 0 */
1675 ptaAddPt(ptad, x, y);
1676 }
1677 numaDestroy(&naerror);
1678 n = ptaGetCount(ptad);
1679 if ((n = ptaGetCount(ptad)) < 3) {
1680 ptaDestroy(&ptad);
1681 return ERROR_INT("less than 3 pts found", __func__, 1);
1682 }
1683
1684 /* Do LSF again */
1685 ret = ptaGetQuadraticLSF(ptad, pa, pb, pc, pnafit);
1686 if (pptad)
1687 *pptad = ptad;
1688 else
1689 ptaDestroy(&ptad);
1690
1691 return ret;
1692 }
1693
1694
1695 /*!
1696 * \brief applyLinearFit()
1697 *
1698 * \param[in] a, b linear fit coefficients
1699 * \param[in] x
1700 * \param[out] py y = a * x + b
1701 * \return 0 if OK, 1 on error
1702 */
1703 l_ok
1704 applyLinearFit(l_float32 a,
1705 l_float32 b,
1706 l_float32 x,
1707 l_float32 *py)
1708 {
1709 if (!py)
1710 return ERROR_INT("&y not defined", __func__, 1);
1711
1712 *py = a * x + b;
1713 return 0;
1714 }
1715
1716
1717 /*!
1718 * \brief applyQuadraticFit()
1719 *
1720 * \param[in] a, b, c quadratic fit coefficients
1721 * \param[in] x
1722 * \param[out] py y = a * x^2 + b * x + c
1723 * \return 0 if OK, 1 on error
1724 */
1725 l_ok
1726 applyQuadraticFit(l_float32 a,
1727 l_float32 b,
1728 l_float32 c,
1729 l_float32 x,
1730 l_float32 *py)
1731 {
1732 if (!py)
1733 return ERROR_INT("&y not defined", __func__, 1);
1734
1735 *py = a * x * x + b * x + c;
1736 return 0;
1737 }
1738
1739
1740 /*!
1741 * \brief applyCubicFit()
1742 *
1743 * \param[in] a, b, c, d cubic fit coefficients
1744 * \param[in] x
1745 * \param[out] py y = a * x^3 + b * x^2 + c * x + d
1746 * \return 0 if OK, 1 on error
1747 */
1748 l_ok
1749 applyCubicFit(l_float32 a,
1750 l_float32 b,
1751 l_float32 c,
1752 l_float32 d,
1753 l_float32 x,
1754 l_float32 *py)
1755 {
1756 if (!py)
1757 return ERROR_INT("&y not defined", __func__, 1);
1758
1759 *py = a * x * x * x + b * x * x + c * x + d;
1760 return 0;
1761 }
1762
1763
1764 /*!
1765 * \brief applyQuarticFit()
1766 *
1767 * \param[in] a, b, c, d, e quartic fit coefficients
1768 * \param[in] x
1769 * \param[out] py y = a * x^4 + b * x^3 + c * x^2 + d * x + e
1770 * \return 0 if OK, 1 on error
1771 */
1772 l_ok
1773 applyQuarticFit(l_float32 a,
1774 l_float32 b,
1775 l_float32 c,
1776 l_float32 d,
1777 l_float32 e,
1778 l_float32 x,
1779 l_float32 *py)
1780 {
1781 l_float32 x2;
1782
1783 if (!py)
1784 return ERROR_INT("&y not defined", __func__, 1);
1785
1786 x2 = x * x;
1787 *py = a * x2 * x2 + b * x2 * x + c * x2 + d * x + e;
1788 return 0;
1789 }
1790
1791
1792 /*---------------------------------------------------------------------*
1793 * Interconversions with Pix *
1794 *---------------------------------------------------------------------*/
1795 /*!
1796 * \brief pixPlotAlongPta()
1797 *
1798 * \param[in] pixs any depth
1799 * \param[in] pta set of points on which to plot
1800 * \param[in] outformat GPLOT_PNG, GPLOT_PS, GPLOT_EPS, GPLOT_LATEX
1801 * \param[in] title [optional] for plot; can be null
1802 * \return 0 if OK, 1 on error
1803 *
1804 * <pre>
1805 * Notes:
1806 * (1) This is a debugging function.
1807 * (2) Removes existing colormaps and clips the pta to the input %pixs.
1808 * (3) If the image is RGB, three separate plots are generated.
1809 * </pre>
1810 */
1811 l_ok
1812 pixPlotAlongPta(PIX *pixs,
1813 PTA *pta,
1814 l_int32 outformat,
1815 const char *title)
1816 {
1817 char buffer[128];
1818 char *rtitle, *gtitle, *btitle;
1819 static l_int32 count = 0; /* require separate temp files for each call */
1820 l_int32 i, x, y, d, w, h, npts, rval, gval, bval;
1821 l_uint32 val;
1822 NUMA *na, *nar, *nag, *nab;
1823 PIX *pixt;
1824
1825 lept_mkdir("lept/plot");
1826
1827 if (!pixs)
1828 return ERROR_INT("pixs not defined", __func__, 1);
1829 if (!pta)
1830 return ERROR_INT("pta not defined", __func__, 1);
1831 if (outformat != GPLOT_PNG && outformat != GPLOT_PS &&
1832 outformat != GPLOT_EPS && outformat != GPLOT_LATEX) {
1833 L_WARNING("outformat invalid; using GPLOT_PNG\n", __func__);
1834 outformat = GPLOT_PNG;
1835 }
1836
1837 pixt = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
1838 d = pixGetDepth(pixt);
1839 w = pixGetWidth(pixt);
1840 h = pixGetHeight(pixt);
1841 npts = ptaGetCount(pta);
1842 if (d == 32) {
1843 nar = numaCreate(npts);
1844 nag = numaCreate(npts);
1845 nab = numaCreate(npts);
1846 for (i = 0; i < npts; i++) {
1847 ptaGetIPt(pta, i, &x, &y);
1848 if (x < 0 || x >= w)
1849 continue;
1850 if (y < 0 || y >= h)
1851 continue;
1852 pixGetPixel(pixt, x, y, &val);
1853 rval = GET_DATA_BYTE(&val, COLOR_RED);
1854 gval = GET_DATA_BYTE(&val, COLOR_GREEN);
1855 bval = GET_DATA_BYTE(&val, COLOR_BLUE);
1856 numaAddNumber(nar, rval);
1857 numaAddNumber(nag, gval);
1858 numaAddNumber(nab, bval);
1859 }
1860
1861 snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
1862 rtitle = stringJoin("Red: ", title);
1863 gplotSimple1(nar, outformat, buffer, rtitle);
1864 snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
1865 gtitle = stringJoin("Green: ", title);
1866 gplotSimple1(nag, outformat, buffer, gtitle);
1867 snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
1868 btitle = stringJoin("Blue: ", title);
1869 gplotSimple1(nab, outformat, buffer, btitle);
1870 numaDestroy(&nar);
1871 numaDestroy(&nag);
1872 numaDestroy(&nab);
1873 LEPT_FREE(rtitle);
1874 LEPT_FREE(gtitle);
1875 LEPT_FREE(btitle);
1876 } else {
1877 na = numaCreate(npts);
1878 for (i = 0; i < npts; i++) {
1879 ptaGetIPt(pta, i, &x, &y);
1880 if (x < 0 || x >= w)
1881 continue;
1882 if (y < 0 || y >= h)
1883 continue;
1884 pixGetPixel(pixt, x, y, &val);
1885 numaAddNumber(na, (l_float32)val);
1886 }
1887
1888 snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
1889 gplotSimple1(na, outformat, buffer, title);
1890 numaDestroy(&na);
1891 }
1892 pixDestroy(&pixt);
1893 return 0;
1894 }
1895
1896
1897 /*!
1898 * \brief ptaGetPixelsFromPix()
1899 *
1900 * \param[in] pixs 1 bpp
1901 * \param[in] box [optional] can be null
1902 * \return pta, or NULL on error
1903 *
1904 * <pre>
1905 * Notes:
1906 * (1) Generates a pta of fg pixels in the pix, within the box.
1907 * If box == NULL, it uses the entire pix.
1908 * </pre>
1909 */
1910 PTA *
1911 ptaGetPixelsFromPix(PIX *pixs,
1912 BOX *box)
1913 {
1914 l_int32 i, j, w, h, wpl, xstart, xend, ystart, yend, bw, bh;
1915 l_uint32 *data, *line;
1916 PTA *pta;
1917
1918 if (!pixs || (pixGetDepth(pixs) != 1))
1919 return (PTA *)ERROR_PTR("pixs undefined or not 1 bpp", __func__, NULL);
1920
1921 pixGetDimensions(pixs, &w, &h, NULL);
1922 data = pixGetData(pixs);
1923 wpl = pixGetWpl(pixs);
1924 xstart = ystart = 0;
1925 xend = w - 1;
1926 yend = h - 1;
1927 if (box) {
1928 boxGetGeometry(box, &xstart, &ystart, &bw, &bh);
1929 xend = xstart + bw - 1;
1930 yend = ystart + bh - 1;
1931 }
1932
1933 if ((pta = ptaCreate(0)) == NULL)
1934 return (PTA *)ERROR_PTR("pta not made", __func__, NULL);
1935 for (i = ystart; i <= yend; i++) {
1936 line = data + i * wpl;
1937 for (j = xstart; j <= xend; j++) {
1938 if (GET_DATA_BIT(line, j))
1939 ptaAddPt(pta, j, i);
1940 }
1941 }
1942
1943 return pta;
1944 }
1945
1946
1947 /*!
1948 * \brief pixGenerateFromPta()
1949 *
1950 * \param[in] pta
1951 * \param[in] w, h of pix
1952 * \return pix 1 bpp, or NULL on error
1953 *
1954 * <pre>
1955 * Notes:
1956 * (1) Points are rounded to nearest ints.
1957 * (2) Any points outside (w,h) are silently discarded.
1958 * (3) Output 1 bpp pix has values 1 for each point in the pta.
1959 * </pre>
1960 */
1961 PIX *
1962 pixGenerateFromPta(PTA *pta,
1963 l_int32 w,
1964 l_int32 h)
1965 {
1966 l_int32 n, i, x, y;
1967 PIX *pix;
1968
1969 if (!pta)
1970 return (PIX *)ERROR_PTR("pta not defined", __func__, NULL);
1971
1972 if ((pix = pixCreate(w, h, 1)) == NULL)
1973 return (PIX *)ERROR_PTR("pix not made", __func__, NULL);
1974 n = ptaGetCount(pta);
1975 for (i = 0; i < n; i++) {
1976 ptaGetIPt(pta, i, &x, &y);
1977 if (x < 0 || x >= w || y < 0 || y >= h)
1978 continue;
1979 pixSetPixel(pix, x, y, 1);
1980 }
1981
1982 return pix;
1983 }
1984
1985
1986 /*!
1987 * \brief ptaGetBoundaryPixels()
1988 *
1989 * \param[in] pixs 1 bpp
1990 * \param[in] type L_BOUNDARY_FG, L_BOUNDARY_BG
1991 * \return pta, or NULL on error
1992 *
1993 * <pre>
1994 * Notes:
1995 * (1) This generates a pta of either fg or bg boundary pixels.
1996 * (2) See also pixGeneratePtaBoundary() for rendering of
1997 * fg boundary pixels.
1998 * </pre>
1999 */
2000 PTA *
2001 ptaGetBoundaryPixels(PIX *pixs,
2002 l_int32 type)
2003 {
2004 PIX *pixt;
2005 PTA *pta;
2006
2007 if (!pixs || (pixGetDepth(pixs) != 1))
2008 return (PTA *)ERROR_PTR("pixs undefined or not 1 bpp", __func__, NULL);
2009 if (type != L_BOUNDARY_FG && type != L_BOUNDARY_BG)
2010 return (PTA *)ERROR_PTR("invalid type", __func__, NULL);
2011
2012 if (type == L_BOUNDARY_FG)
2013 pixt = pixMorphSequence(pixs, "e3.3", 0);
2014 else
2015 pixt = pixMorphSequence(pixs, "d3.3", 0);
2016 pixXor(pixt, pixt, pixs);
2017 pta = ptaGetPixelsFromPix(pixt, NULL);
2018
2019 pixDestroy(&pixt);
2020 return pta;
2021 }
2022
2023
2024 /*!
2025 * \brief ptaaGetBoundaryPixels()
2026 *
2027 * \param[in] pixs 1 bpp
2028 * \param[in] type L_BOUNDARY_FG, L_BOUNDARY_BG
2029 * \param[in] connectivity 4 or 8
2030 * \param[out] pboxa [optional] bounding boxes of the c.c.
2031 * \param[out] ppixa [optional] pixa of the c.c.
2032 * \return ptaa, or NULL on error
2033 *
2034 * <pre>
2035 * Notes:
2036 * (1) This generates a ptaa of either fg or bg boundary pixels,
2037 * where each pta has the boundary pixels for a connected
2038 * component.
2039 * (2) We can't simply find all the boundary pixels and then select
2040 * those within the bounding box of each component, because
2041 * bounding boxes can overlap. It is necessary to extract and
2042 * dilate or erode each component separately. Note also that
2043 * special handling is required for bg pixels when the
2044 * component touches the pix boundary.
2045 * </pre>
2046 */
2047 PTAA *
2048 ptaaGetBoundaryPixels(PIX *pixs,
2049 l_int32 type,
2050 l_int32 connectivity,
2051 BOXA **pboxa,
2052 PIXA **ppixa)
2053 {
2054 l_int32 i, n, w, h, x, y, bw, bh, left, right, top, bot;
2055 BOXA *boxa;
2056 PIX *pixt1, *pixt2;
2057 PIXA *pixa;
2058 PTA *pta1, *pta2;
2059 PTAA *ptaa;
2060
2061 if (pboxa) *pboxa = NULL;
2062 if (ppixa) *ppixa = NULL;
2063 if (!pixs || (pixGetDepth(pixs) != 1))
2064 return (PTAA *)ERROR_PTR("pixs undefined or not 1 bpp", __func__, NULL);
2065 if (type != L_BOUNDARY_FG && type != L_BOUNDARY_BG)
2066 return (PTAA *)ERROR_PTR("invalid type", __func__, NULL);
2067 if (connectivity != 4 && connectivity != 8)
2068 return (PTAA *)ERROR_PTR("connectivity not 4 or 8", __func__, NULL);
2069
2070 pixGetDimensions(pixs, &w, &h, NULL);
2071 boxa = pixConnComp(pixs, &pixa, connectivity);
2072 n = boxaGetCount(boxa);
2073 ptaa = ptaaCreate(0);
2074 for (i = 0; i < n; i++) {
2075 pixt1 = pixaGetPix(pixa, i, L_CLONE);
2076 boxaGetBoxGeometry(boxa, i, &x, &y, &bw, &bh);
2077 left = right = top = bot = 0;
2078 if (type == L_BOUNDARY_BG) {
2079 if (x > 0) left = 1;
2080 if (y > 0) top = 1;
2081 if (x + bw < w) right = 1;
2082 if (y + bh < h) bot = 1;
2083 pixt2 = pixAddBorderGeneral(pixt1, left, right, top, bot, 0);
2084 } else {
2085 pixt2 = pixClone(pixt1);
2086 }
2087 pta1 = ptaGetBoundaryPixels(pixt2, type);
2088 pta2 = ptaTransform(pta1, x - left, y - top, 1.0, 1.0);
2089 ptaaAddPta(ptaa, pta2, L_INSERT);
2090 ptaDestroy(&pta1);
2091 pixDestroy(&pixt1);
2092 pixDestroy(&pixt2);
2093 }
2094
2095 if (pboxa)
2096 *pboxa = boxa;
2097 else
2098 boxaDestroy(&boxa);
2099 if (ppixa)
2100 *ppixa = pixa;
2101 else
2102 pixaDestroy(&pixa);
2103 return ptaa;
2104 }
2105
2106
2107 /*!
2108 * \brief ptaaIndexLabeledPixels()
2109 *
2110 * \param[in] pixs 32 bpp, of indices of c.c.
2111 * \param[out] pncc [optional] number of connected components
2112 * \return ptaa, or NULL on error
2113 *
2114 * <pre>
2115 * Notes:
2116 * (1) The pixel values in %pixs are the index of the connected component
2117 * to which the pixel belongs; %pixs is typically generated from
2118 * a 1 bpp pix by pixConnCompTransform(). Background pixels in
2119 * the generating 1 bpp pix are represented in %pixs by 0.
2120 * We do not check that the pixel values are correctly labelled.
2121 * (2) Each pta in the returned ptaa gives the pixel locations
2122 * corresponding to a connected component, with the label of each
2123 * given by the index of the pta into the ptaa.
2124 * (3) Initialize with the first pta in ptaa being empty and
2125 * representing the background value (index 0) in the pix.
2126 * </pre>
2127 */
2128 PTAA *
2129 ptaaIndexLabeledPixels(PIX *pixs,
2130 l_int32 *pncc)
2131 {
2132 l_int32 wpl, index, i, j, w, h;
2133 l_uint32 maxval;
2134 l_uint32 *data, *line;
2135 PTA *pta;
2136 PTAA *ptaa;
2137
2138 if (pncc) *pncc = 0;
2139 if (!pixs || (pixGetDepth(pixs) != 32))
2140 return (PTAA *)ERROR_PTR("pixs undef or not 32 bpp", __func__, NULL);
2141
2142 /* The number of c.c. is the maximum pixel value. Use this to
2143 * initialize ptaa with sufficient pta arrays */
2144 pixGetMaxValueInRect(pixs, NULL, &maxval, NULL, NULL);
2145 if (pncc) *pncc = maxval;
2146 pta = ptaCreate(1);
2147 ptaa = ptaaCreate(maxval + 1);
2148 ptaaInitFull(ptaa, pta);
2149 ptaDestroy(&pta);
2150
2151 /* Sweep over %pixs, saving the pixel coordinates of each pixel
2152 * with nonzero value in the appropriate pta, indexed by that value. */
2153 pixGetDimensions(pixs, &w, &h, NULL);
2154 data = pixGetData(pixs);
2155 wpl = pixGetWpl(pixs);
2156 for (i = 0; i < h; i++) {
2157 line = data + wpl * i;
2158 for (j = 0; j < w; j++) {
2159 index = line[j];
2160 if (index > 0)
2161 ptaaAddPt(ptaa, index, j, i);
2162 }
2163 }
2164
2165 return ptaa;
2166 }
2167
2168
2169 /*!
2170 * \brief ptaGetNeighborPixLocs()
2171 *
2172 * \param[in] pixs any depth
2173 * \param[in] x, y pixel from which we search for nearest neighbors
2174 * \param[in] conn 4 or 8 connectivity
2175 * \return pta, or NULL on error
2176 *
2177 * <pre>
2178 * Notes:
2179 * (1) Generates a pta of all valid neighbor pixel locations,
2180 * or NULL on error.
2181 * </pre>
2182 */
2183 PTA *
2184 ptaGetNeighborPixLocs(PIX *pixs,
2185 l_int32 x,
2186 l_int32 y,
2187 l_int32 conn)
2188 {
2189 l_int32 w, h;
2190 PTA *pta;
2191
2192 if (!pixs)
2193 return (PTA *)ERROR_PTR("pixs not defined", __func__, NULL);
2194 pixGetDimensions(pixs, &w, &h, NULL);
2195 if (x < 0 || x >= w || y < 0 || y >= h)
2196 return (PTA *)ERROR_PTR("(x,y) not in pixs", __func__, NULL);
2197 if (conn != 4 && conn != 8)
2198 return (PTA *)ERROR_PTR("conn not 4 or 8", __func__, NULL);
2199
2200 pta = ptaCreate(conn);
2201 if (x > 0)
2202 ptaAddPt(pta, x - 1, y);
2203 if (x < w - 1)
2204 ptaAddPt(pta, x + 1, y);
2205 if (y > 0)
2206 ptaAddPt(pta, x, y - 1);
2207 if (y < h - 1)
2208 ptaAddPt(pta, x, y + 1);
2209 if (conn == 8) {
2210 if (x > 0) {
2211 if (y > 0)
2212 ptaAddPt(pta, x - 1, y - 1);
2213 if (y < h - 1)
2214 ptaAddPt(pta, x - 1, y + 1);
2215 }
2216 if (x < w - 1) {
2217 if (y > 0)
2218 ptaAddPt(pta, x + 1, y - 1);
2219 if (y < h - 1)
2220 ptaAddPt(pta, x + 1, y + 1);
2221 }
2222 }
2223
2224 return pta;
2225 }
2226
2227
2228 /*---------------------------------------------------------------------*
2229 * Interconversion with Numa *
2230 *---------------------------------------------------------------------*/
2231 /*!
2232 * \brief numaConvertToPta1()
2233 *
2234 * \param[in] na numa with implicit y(x)
2235 * \return pta if OK; null on error
2236 */
2237 PTA *
2238 numaConvertToPta1(NUMA *na)
2239 {
2240 l_int32 i, n;
2241 l_float32 startx, delx, val;
2242 PTA *pta;
2243
2244 if (!na)
2245 return (PTA *)ERROR_PTR("na not defined", __func__, NULL);
2246
2247 n = numaGetCount(na);
2248 pta = ptaCreate(n);
2249 numaGetParameters(na, &startx, &delx);
2250 for (i = 0; i < n; i++) {
2251 numaGetFValue(na, i, &val);
2252 ptaAddPt(pta, startx + i * delx, val);
2253 }
2254 return pta;
2255 }
2256
2257
2258 /*!
2259 * \brief numaConvertToPta2()
2260 *
2261 * \param[in] nax
2262 * \param[in] nay
2263 * \return pta if OK; null on error
2264 */
2265 PTA *
2266 numaConvertToPta2(NUMA *nax,
2267 NUMA *nay)
2268 {
2269 l_int32 i, n, nx, ny;
2270 l_float32 valx, valy;
2271 PTA *pta;
2272
2273 if (!nax || !nay)
2274 return (PTA *)ERROR_PTR("nax and nay not both defined", __func__, NULL);
2275
2276 nx = numaGetCount(nax);
2277 ny = numaGetCount(nay);
2278 n = L_MIN(nx, ny);
2279 if (nx != ny)
2280 L_WARNING("nx = %d does not equal ny = %d\n", __func__, nx, ny);
2281 pta = ptaCreate(n);
2282 for (i = 0; i < n; i++) {
2283 numaGetFValue(nax, i, &valx);
2284 numaGetFValue(nay, i, &valy);
2285 ptaAddPt(pta, valx, valy);
2286 }
2287 return pta;
2288 }
2289
2290
2291 /*!
2292 * \brief ptaConvertToNuma()
2293 *
2294 * \param[in] pta
2295 * \param[out] pnax addr of nax
2296 * \param[out] pnay addr of nay
2297 * \return 0 if OK, 1 on error
2298 */
2299 l_ok
2300 ptaConvertToNuma(PTA *pta,
2301 NUMA **pnax,
2302 NUMA **pnay)
2303 {
2304 l_int32 i, n;
2305 l_float32 valx, valy;
2306
2307 if (pnax) *pnax = NULL;
2308 if (pnay) *pnay = NULL;
2309 if (!pnax || !pnay)
2310 return ERROR_INT("&nax and &nay not both defined", __func__, 1);
2311 if (!pta)
2312 return ERROR_INT("pta not defined", __func__, 1);
2313
2314 n = ptaGetCount(pta);
2315 *pnax = numaCreate(n);
2316 *pnay = numaCreate(n);
2317 for (i = 0; i < n; i++) {
2318 ptaGetPt(pta, i, &valx, &valy);
2319 numaAddNumber(*pnax, valx);
2320 numaAddNumber(*pnay, valy);
2321 }
2322 return 0;
2323 }
2324
2325
2326 /*---------------------------------------------------------------------*
2327 * Display Pta and Ptaa *
2328 *---------------------------------------------------------------------*/
2329 /*!
2330 * \brief pixDisplayPta()
2331 *
2332 * \param[in] pixd can be same as pixs or NULL; 32 bpp if in-place
2333 * \param[in] pixs 1, 2, 4, 8, 16 or 32 bpp
2334 * \param[in] pta of path to be plotted
2335 * \return pixd 32 bpp RGB version of pixs, with path in green.
2336 *
2337 * <pre>
2338 * Notes:
2339 * (1) To write on an existing pixs, pixs must be 32 bpp and
2340 * call with pixd == pixs:
2341 * pixDisplayPta(pixs, pixs, pta);
2342 * To write to a new pix, use pixd == NULL and call:
2343 * pixd = pixDisplayPta(NULL, pixs, pta);
2344 * (2) On error, returns pixd to avoid losing pixs if called as
2345 * pixs = pixDisplayPta(pixs, pixs, pta);
2346 * </pre>
2347 */
2348 PIX *
2349 pixDisplayPta(PIX *pixd,
2350 PIX *pixs,
2351 PTA *pta)
2352 {
2353 l_int32 i, n, w, h, x, y;
2354 l_uint32 rpixel, gpixel, bpixel;
2355
2356 if (!pixs)
2357 return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
2358 if (!pta)
2359 return (PIX *)ERROR_PTR("pta not defined", __func__, pixd);
2360 if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
2361 return (PIX *)ERROR_PTR("invalid pixd", __func__, pixd);
2362
2363 if (!pixd)
2364 pixd = pixConvertTo32(pixs);
2365 pixGetDimensions(pixd, &w, &h, NULL);
2366 composeRGBPixel(255, 0, 0, &rpixel); /* start point */
2367 composeRGBPixel(0, 255, 0, &gpixel);
2368 composeRGBPixel(0, 0, 255, &bpixel); /* end point */
2369
2370 n = ptaGetCount(pta);
2371 for (i = 0; i < n; i++) {
2372 ptaGetIPt(pta, i, &x, &y);
2373 if (x < 0 || x >= w || y < 0 || y >= h)
2374 continue;
2375 if (i == 0)
2376 pixSetPixel(pixd, x, y, rpixel);
2377 else if (i < n - 1)
2378 pixSetPixel(pixd, x, y, gpixel);
2379 else
2380 pixSetPixel(pixd, x, y, bpixel);
2381 }
2382
2383 return pixd;
2384 }
2385
2386
2387 /*!
2388 * \brief pixDisplayPtaaPattern()
2389 *
2390 * \param[in] pixd 32 bpp
2391 * \param[in] pixs 1, 2, 4, 8, 16 or 32 bpp; 32 bpp if in place
2392 * \param[in] ptaa giving locations at which the pattern is displayed
2393 * \param[in] pixp 1 bpp pattern to be placed such that its reference
2394 * point co-locates with each point in pta
2395 * \param[in] cx, cy reference point in pattern
2396 * \return pixd 32 bpp RGB version of pixs.
2397 *
2398 * <pre>
2399 * Notes:
2400 * (1) To write on an existing pixs, pixs must be 32 bpp and
2401 * call with pixd == pixs:
2402 * pixDisplayPtaPattern(pixs, pixs, pta, ...);
2403 * To write to a new pix, use pixd == NULL and call:
2404 * pixd = pixDisplayPtaPattern(NULL, pixs, pta, ...);
2405 * (2) Puts a random color on each pattern associated with a pta.
2406 * (3) On error, returns pixd to avoid losing pixs if called as
2407 * pixs = pixDisplayPtaPattern(pixs, pixs, pta, ...);
2408 * (4) A typical pattern to be used is a circle, generated with
2409 * generatePtaFilledCircle()
2410 * </pre>
2411 */
2412 PIX *
2413 pixDisplayPtaaPattern(PIX *pixd,
2414 PIX *pixs,
2415 PTAA *ptaa,
2416 PIX *pixp,
2417 l_int32 cx,
2418 l_int32 cy)
2419 {
2420 l_int32 i, n;
2421 l_uint32 color;
2422 PIXCMAP *cmap;
2423 PTA *pta;
2424
2425 if (!pixs)
2426 return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
2427 if (!ptaa)
2428 return (PIX *)ERROR_PTR("ptaa not defined", __func__, pixd);
2429 if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
2430 return (PIX *)ERROR_PTR("invalid pixd", __func__, pixd);
2431 if (!pixp)
2432 return (PIX *)ERROR_PTR("pixp not defined", __func__, pixd);
2433
2434 if (!pixd)
2435 pixd = pixConvertTo32(pixs);
2436
2437 /* Use 256 random colors */
2438 cmap = pixcmapCreateRandom(8, 0, 0);
2439 n = ptaaGetCount(ptaa);
2440 for (i = 0; i < n; i++) {
2441 pixcmapGetColor32(cmap, i % 256, &color);
2442 pta = ptaaGetPta(ptaa, i, L_CLONE);
2443 pixDisplayPtaPattern(pixd, pixd, pta, pixp, cx, cy, color);
2444 ptaDestroy(&pta);
2445 }
2446
2447 pixcmapDestroy(&cmap);
2448 return pixd;
2449 }
2450
2451
2452 /*!
2453 * \brief pixDisplayPtaPattern()
2454 *
2455 * \param[in] pixd can be same as pixs or NULL; 32 bpp if in-place
2456 * \param[in] pixs 1, 2, 4, 8, 16 or 32 bpp
2457 * \param[in] pta giving locations at which the pattern is displayed
2458 * \param[in] pixp 1 bpp pattern to be placed such that its reference
2459 * point co-locates with each point in pta
2460 * \param[in] cx, cy reference point in pattern
2461 * \param[in] color in 0xrrggbb00 format
2462 * \return pixd 32 bpp RGB version of pixs.
2463 *
2464 * <pre>
2465 * Notes:
2466 * (1) To write on an existing pixs, pixs must be 32 bpp and
2467 * call with pixd == pixs:
2468 * pixDisplayPtaPattern(pixs, pixs, pta, ...);
2469 * To write to a new pix, use pixd == NULL and call:
2470 * pixd = pixDisplayPtaPattern(NULL, pixs, pta, ...);
2471 * (2) On error, returns pixd to avoid losing pixs if called as
2472 * pixs = pixDisplayPtaPattern(pixs, pixs, pta, ...);
2473 * (3) A typical pattern to be used is a circle, generated with
2474 * generatePtaFilledCircle()
2475 * </pre>
2476 */
2477 PIX *
2478 pixDisplayPtaPattern(PIX *pixd,
2479 PIX *pixs,
2480 PTA *pta,
2481 PIX *pixp,
2482 l_int32 cx,
2483 l_int32 cy,
2484 l_uint32 color)
2485 {
2486 l_int32 i, n, w, h, x, y;
2487 PTA *ptat;
2488
2489 if (!pixs)
2490 return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
2491 if (!pta)
2492 return (PIX *)ERROR_PTR("pta not defined", __func__, pixd);
2493 if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
2494 return (PIX *)ERROR_PTR("invalid pixd", __func__, pixd);
2495 if (!pixp)
2496 return (PIX *)ERROR_PTR("pixp not defined", __func__, pixd);
2497
2498 if (!pixd)
2499 pixd = pixConvertTo32(pixs);
2500 pixGetDimensions(pixs, &w, &h, NULL);
2501 ptat = ptaReplicatePattern(pta, pixp, NULL, cx, cy, w, h);
2502
2503 n = ptaGetCount(ptat);
2504 for (i = 0; i < n; i++) {
2505 ptaGetIPt(ptat, i, &x, &y);
2506 if (x < 0 || x >= w || y < 0 || y >= h)
2507 continue;
2508 pixSetPixel(pixd, x, y, color);
2509 }
2510
2511 ptaDestroy(&ptat);
2512 return pixd;
2513 }
2514
2515
2516 /*!
2517 * \brief ptaReplicatePattern()
2518 *
2519 * \param[in] ptas "sparse" input pta
2520 * \param[in] pixp [optional] 1 bpp pattern, to be replicated
2521 * in output pta
2522 * \param[in] ptap [optional] set of pts, to be replicated in output pta
2523 * \param[in] cx, cy reference point in pattern
2524 * \param[in] w, h clipping sizes for output pta
2525 * \return ptad with all points of replicated pattern, or NULL on error
2526 *
2527 * <pre>
2528 * Notes:
2529 * (1) You can use either the image %pixp or the set of pts %ptap.
2530 * (2) The pattern is placed with its reference point at each point
2531 * in ptas, and all the fg pixels are colleced into ptad.
2532 * For %pixp, this is equivalent to blitting pixp at each point
2533 * in ptas, and then converting the resulting pix to a pta.
2534 * </pre>
2535 */
2536 PTA *
2537 ptaReplicatePattern(PTA *ptas,
2538 PIX *pixp,
2539 PTA *ptap,
2540 l_int32 cx,
2541 l_int32 cy,
2542 l_int32 w,
2543 l_int32 h)
2544 {
2545 l_int32 i, j, n, np, x, y, xp, yp, xf, yf;
2546 PTA *ptat, *ptad;
2547
2548 if (!ptas)
2549 return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
2550 if (!pixp && !ptap)
2551 return (PTA *)ERROR_PTR("no pattern is defined", __func__, NULL);
2552 if (pixp && ptap)
2553 L_WARNING("pixp and ptap defined; using ptap\n", __func__);
2554
2555 n = ptaGetCount(ptas);
2556 ptad = ptaCreate(n);
2557 if (ptap)
2558 ptat = ptaClone(ptap);
2559 else
2560 ptat = ptaGetPixelsFromPix(pixp, NULL);
2561 np = ptaGetCount(ptat);
2562 for (i = 0; i < n; i++) {
2563 ptaGetIPt(ptas, i, &x, &y);
2564 for (j = 0; j < np; j++) {
2565 ptaGetIPt(ptat, j, &xp, &yp);
2566 xf = x - cx + xp;
2567 yf = y - cy + yp;
2568 if (xf >= 0 && xf < w && yf >= 0 && yf < h)
2569 ptaAddPt(ptad, xf, yf);
2570 }
2571 }
2572
2573 ptaDestroy(&ptat);
2574 return ptad;
2575 }
2576
2577
2578 /*!
2579 * \brief pixDisplayPtaa()
2580 *
2581 * \param[in] pixs 1, 2, 4, 8, 16 or 32 bpp
2582 * \param[in] ptaa array of paths to be plotted
2583 * \return pixd 32 bpp RGB version of pixs, with paths plotted
2584 * in different colors, or NULL on error
2585 */
2586 PIX *
2587 pixDisplayPtaa(PIX *pixs,
2588 PTAA *ptaa)
2589 {
2590 l_int32 i, j, w, h, npta, npt, x, y, rv, gv, bv;
2591 l_uint32 *pixela;
2592 NUMA *na1, *na2, *na3;
2593 PIX *pixd;
2594 PTA *pta;
2595
2596 if (!pixs)
2597 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
2598 if (!ptaa)
2599 return (PIX *)ERROR_PTR("ptaa not defined", __func__, NULL);
2600 npta = ptaaGetCount(ptaa);
2601 if (npta == 0)
2602 return (PIX *)ERROR_PTR("no pta", __func__, NULL);
2603
2604 if ((pixd = pixConvertTo32(pixs)) == NULL)
2605 return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
2606 pixGetDimensions(pixd, &w, &h, NULL);
2607
2608 /* Make a colormap for the paths */
2609 if ((pixela = (l_uint32 *)LEPT_CALLOC(npta, sizeof(l_uint32))) == NULL) {
2610 pixDestroy(&pixd);
2611 return (PIX *)ERROR_PTR("calloc fail for pixela", __func__, NULL);
2612 }
2613 na1 = numaPseudorandomSequence(256, 14657);
2614 na2 = numaPseudorandomSequence(256, 34631);
2615 na3 = numaPseudorandomSequence(256, 54617);
2616 for (i = 0; i < npta; i++) {
2617 numaGetIValue(na1, i % 256, &rv);
2618 numaGetIValue(na2, i % 256, &gv);
2619 numaGetIValue(na3, i % 256, &bv);
2620 composeRGBPixel(rv, gv, bv, &pixela[i]);
2621 }
2622 numaDestroy(&na1);
2623 numaDestroy(&na2);
2624 numaDestroy(&na3);
2625
2626 for (i = 0; i < npta; i++) {
2627 pta = ptaaGetPta(ptaa, i, L_CLONE);
2628 npt = ptaGetCount(pta);
2629 for (j = 0; j < npt; j++) {
2630 ptaGetIPt(pta, j, &x, &y);
2631 if (x < 0 || x >= w || y < 0 || y >= h)
2632 continue;
2633 pixSetPixel(pixd, x, y, pixela[i]);
2634 }
2635 ptaDestroy(&pta);
2636 }
2637
2638 LEPT_FREE(pixela);
2639 return pixd;
2640 }