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