Coverage Report

Created: 2025-06-13 07:15

/src/leptonica/src/ptafunc1.c
Line
Count
Source (jump to first uncovered line)
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
0
{
128
0
l_int32    n, i;
129
0
l_float32  x, y;
130
0
PTA       *ptad;
131
132
0
    if (!ptas)
133
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
134
0
    if (subfactor < 1)
135
0
        return (PTA *)ERROR_PTR("subfactor < 1", __func__, NULL);
136
137
0
    ptad = ptaCreate(0);
138
0
    n = ptaGetCount(ptas);
139
0
    for (i = 0; i < n; i++) {
140
0
        if (i % subfactor != 0) continue;
141
0
        ptaGetPt(ptas, i, &x, &y);
142
0
        ptaAddPt(ptad, x, y);
143
0
    }
144
145
0
    return ptad;
146
0
}
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
0
{
171
0
l_int32  n, i, x, y;
172
173
0
    if (!ptad)
174
0
        return ERROR_INT("ptad not defined", __func__, 1);
175
0
    if (!ptas)
176
0
        return 0;
177
178
0
    if (istart < 0)
179
0
        istart = 0;
180
0
    n = ptaGetCount(ptas);
181
0
    if (iend < 0 || iend >= n)
182
0
        iend = n - 1;
183
0
    if (istart > iend)
184
0
        return ERROR_INT("istart > iend; no pts", __func__, 1);
185
186
0
    for (i = istart; i <= iend; i++) {
187
0
        ptaGetIPt(ptas, i, &x, &y);
188
0
        if (ptaAddPt(ptad, x, y) == 1) {
189
0
            L_ERROR("failed to add pt at i = %d\n", __func__, i);
190
0
            return 1;
191
0
        }
192
0
    }
193
0
    return 0;
194
0
}
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
0
{
219
0
l_int32  n, i;
220
0
PTA     *pta;
221
222
0
    if (!ptaad)
223
0
        return ERROR_INT("ptaad not defined", __func__, 1);
224
0
    if (!ptaas)
225
0
        return 0;
226
227
0
    if (istart < 0)
228
0
        istart = 0;
229
0
    n = ptaaGetCount(ptaas);
230
0
    if (iend < 0 || iend >= n)
231
0
        iend = n - 1;
232
0
    if (istart > iend)
233
0
        return ERROR_INT("istart > iend; no pts", __func__, 1);
234
235
0
    for (i = istart; i <= iend; i++) {
236
0
        pta = ptaaGetPta(ptaas, i, L_CLONE);
237
0
        ptaaAddPta(ptaad, pta, L_INSERT);
238
0
    }
239
240
0
    return 0;
241
0
}
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
0
{
255
0
l_int32    n, i, ix, iy;
256
0
l_float32  x, y;
257
0
PTA       *ptad;
258
259
0
    if (!ptas)
260
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
261
262
0
    n = ptaGetCount(ptas);
263
0
    if ((ptad = ptaCreate(n)) == NULL)
264
0
        return (PTA *)ERROR_PTR("ptad not made", __func__, NULL);
265
0
    for (i = n - 1; i >= 0; i--) {
266
0
        if (type == 0) {
267
0
            ptaGetPt(ptas, i, &x, &y);
268
0
            ptaAddPt(ptad, x, y);
269
0
        } else {  /* type == 1 */
270
0
            ptaGetIPt(ptas, i, &ix, &iy);
271
0
            ptaAddPt(ptad, ix, iy);
272
0
        }
273
0
    }
274
275
0
    return ptad;
276
0
}
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
0
{
288
0
l_int32    n, i;
289
0
l_float32  x, y;
290
0
PTA       *ptad;
291
292
0
    if (!ptas)
293
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
294
295
0
    n = ptaGetCount(ptas);
296
0
    if ((ptad = ptaCreate(n)) == NULL)
297
0
        return (PTA *)ERROR_PTR("ptad not made", __func__, NULL);
298
0
    for (i = 0; i < n; i++) {
299
0
        ptaGetPt(ptas, i, &x, &y);
300
0
        ptaAddPt(ptad, y, x);
301
0
    }
302
303
0
    return ptad;
304
0
}
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
0
{
328
0
l_int32  n, i, x, y, j, index, state;
329
0
l_int32  x1, y1, x2, y2;
330
0
PTA     *ptad;
331
332
0
    if (!ptas)
333
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
334
335
0
    n = ptaGetCount(ptas);
336
337
        /* Verify input data */
338
0
    ptaGetIPt(ptas, 0, &x1, &y1);
339
0
    ptaGetIPt(ptas, n - 1, &x2, &y2);
340
0
    if (x1 != x2 || y1 != y2)
341
0
        return (PTA *)ERROR_PTR("start and end pts not same", __func__, NULL);
342
0
    state = L_NOT_FOUND;
343
0
    for (i = 0; i < n; i++) {
344
0
        ptaGetIPt(ptas, i, &x, &y);
345
0
        if (x == xs && y == ys) {
346
0
            state = L_FOUND;
347
0
            break;
348
0
        }
349
0
    }
350
0
    if (state == L_NOT_FOUND)
351
0
        return (PTA *)ERROR_PTR("start pt not in ptas", __func__, NULL);
352
353
0
    if ((ptad = ptaCreate(n)) == NULL)
354
0
        return (PTA *)ERROR_PTR("ptad not made", __func__, NULL);
355
0
    for (j = 0; j < n - 1; j++) {
356
0
        if (i + j < n - 1)
357
0
            index = i + j;
358
0
        else
359
0
            index = (i + j + 1) % n;
360
0
        ptaGetIPt(ptas, index, &x, &y);
361
0
        ptaAddPt(ptad, x, y);
362
0
    }
363
0
    ptaAddPt(ptad, xs, ys);
364
365
0
    return ptad;
366
0
}
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
0
{
382
0
l_int32    n, npt, i;
383
0
l_float32  x, y;
384
0
PTA       *ptad;
385
386
0
    if (!ptas)
387
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
388
0
    if ((n = ptaGetCount(ptas)) == 0) {
389
0
        L_WARNING("ptas is empty\n", __func__);
390
0
        return ptaCopy(ptas);
391
0
    }
392
0
    first = L_MAX(0, first);
393
0
    if (last < 0) last = n - 1;
394
0
    if (first >= n)
395
0
        return (PTA *)ERROR_PTR("invalid first", __func__, NULL);
396
0
    if (last >= n) {
397
0
        L_WARNING("last = %d is beyond max index = %d; adjusting\n",
398
0
                  __func__, last, n - 1);
399
0
        last = n - 1;
400
0
    }
401
0
    if (first > last)
402
0
        return (PTA *)ERROR_PTR("first > last", __func__, NULL);
403
404
0
    npt = last - first + 1;
405
0
    ptad = ptaCreate(npt);
406
0
    for (i = first; i <= last; i++) {
407
0
        ptaGetPt(ptas, i, &x, &y);
408
0
        ptaAddPt(ptad, x, y);
409
0
    }
410
0
    return ptad;
411
0
}
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
0
{
433
0
l_int32  n, i, x, y, minx, maxx, miny, maxy;
434
435
0
    if (!pta)
436
0
        return (BOX *)ERROR_PTR("pta not defined", __func__, NULL);
437
438
0
    minx = 10000000;
439
0
    miny = 10000000;
440
0
    maxx = -10000000;
441
0
    maxy = -10000000;
442
0
    n = ptaGetCount(pta);
443
0
    for (i = 0; i < n; i++) {
444
0
        ptaGetIPt(pta, i, &x, &y);
445
0
        if (x < minx) minx = x;
446
0
        if (x > maxx) maxx = x;
447
0
        if (y < miny) miny = y;
448
0
        if (y > maxy) maxy = y;
449
0
    }
450
451
0
    return boxCreate(minx, miny, maxx - minx + 1, maxy - miny + 1);
452
0
}
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
0
{
479
0
l_int32    n, i;
480
0
l_float32  x, y, minx, maxx, miny, maxy;
481
482
0
    if (!pminx && !pmaxx && !pminy && !pmaxy)
483
0
        return ERROR_INT("no output requested", __func__, 1);
484
0
    if (pminx) *pminx = 0;
485
0
    if (pmaxx) *pmaxx = 0;
486
0
    if (pminy) *pminy = 0;
487
0
    if (pmaxy) *pmaxy = 0;
488
0
    if (!pta)
489
0
        return ERROR_INT("pta not defined", __func__, 1);
490
0
    if ((n = ptaGetCount(pta)) == 0)
491
0
        return ERROR_INT("no points in pta", __func__, 1);
492
493
0
    ptaGetPt(pta, 0, &x, &y);
494
0
    minx = x;
495
0
    maxx = x;
496
0
    miny = y;
497
0
    maxy = y;
498
0
    for (i = 1; i < n; i++) {
499
0
        ptaGetPt(pta, i, &x, &y);
500
0
        if (x < minx) minx = x;
501
0
        if (x > maxx) maxx = x;
502
0
        if (y < miny) miny = y;
503
0
        if (y > maxy) maxy = y;
504
0
    }
505
0
    if (pminx) *pminx = minx;
506
0
    if (pmaxx) *pmaxx = maxx;
507
0
    if (pminy) *pminy = miny;
508
0
    if (pmaxy) *pmaxy = maxy;
509
0
    return 0;
510
0
}
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
0
{
524
0
PTA       *ptad;
525
0
l_int32    n, i, contains;
526
0
l_float32  x, y;
527
528
0
    if (!ptas)
529
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
530
0
    if (!box)
531
0
        return (PTA *)ERROR_PTR("box not defined", __func__, NULL);
532
533
0
    n = ptaGetCount(ptas);
534
0
    ptad = ptaCreate(0);
535
0
    for (i = 0; i < n; i++) {
536
0
        ptaGetPt(ptas, i, &x, &y);
537
0
        boxContainsPt(box, x, y, &contains);
538
0
        if (contains)
539
0
            ptaAddPt(ptad, x, y);
540
0
    }
541
542
0
    return ptad;
543
0
}
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
0
{
561
0
l_int32    i, j, x, y, w, h, wpl, mindim, found;
562
0
l_uint32  *data, *line;
563
0
PTA       *pta;
564
565
0
    if (!pixs)
566
0
        return (PTA *)ERROR_PTR("pixs not defined", __func__, NULL);
567
0
    if (pixGetDepth(pixs) != 1)
568
0
        return (PTA *)ERROR_PTR("pixs not 1 bpp", __func__, NULL);
569
570
0
    w = pixGetWidth(pixs);
571
0
    h = pixGetHeight(pixs);
572
0
    mindim = L_MIN(w, h);
573
0
    data = pixGetData(pixs);
574
0
    wpl = pixGetWpl(pixs);
575
576
0
    if ((pta = ptaCreate(4)) == NULL)
577
0
        return (PTA *)ERROR_PTR("pta not made", __func__, NULL);
578
579
0
    for (found = FALSE, i = 0; i < mindim; i++) {
580
0
        for (j = 0; j <= i; j++) {
581
0
            y = i - j;
582
0
            line = data + y * wpl;
583
0
            if (GET_DATA_BIT(line, j)) {
584
0
                ptaAddPt(pta, j, y);
585
0
                found = TRUE;
586
0
                break;
587
0
            }
588
0
        }
589
0
        if (found == TRUE)
590
0
            break;
591
0
    }
592
593
0
    for (found = FALSE, i = 0; i < mindim; i++) {
594
0
        for (j = 0; j <= i; j++) {
595
0
            y = i - j;
596
0
            line = data + y * wpl;
597
0
            x = w - 1 - j;
598
0
            if (GET_DATA_BIT(line, x)) {
599
0
                ptaAddPt(pta, x, y);
600
0
                found = TRUE;
601
0
                break;
602
0
            }
603
0
        }
604
0
        if (found == TRUE)
605
0
            break;
606
0
    }
607
608
0
    for (found = FALSE, i = 0; i < mindim; i++) {
609
0
        for (j = 0; j <= i; j++) {
610
0
            y = h - 1 - i + j;
611
0
            line = data + y * wpl;
612
0
            if (GET_DATA_BIT(line, j)) {
613
0
                ptaAddPt(pta, j, y);
614
0
                found = TRUE;
615
0
                break;
616
0
            }
617
0
        }
618
0
        if (found == TRUE)
619
0
            break;
620
0
    }
621
622
0
    for (found = FALSE, i = 0; i < mindim; i++) {
623
0
        for (j = 0; j <= i; j++) {
624
0
            y = h - 1 - i + j;
625
0
            line = data + y * wpl;
626
0
            x = w - 1 - j;
627
0
            if (GET_DATA_BIT(line, x)) {
628
0
                ptaAddPt(pta, x, y);
629
0
                found = TRUE;
630
0
                break;
631
0
            }
632
0
        }
633
0
        if (found == TRUE)
634
0
            break;
635
0
    }
636
637
0
    return pta;
638
0
}
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
0
{
653
0
l_int32  i, n, ix, iy;
654
655
0
    if (!pta)
656
0
        return ERROR_INT("pta not defined", __func__, 0);
657
658
0
    n = ptaGetCount(pta);
659
0
    for (i = 0; i < n; i++) {
660
0
        ptaGetIPt(pta, i, &ix, &iy);
661
0
        if (x == ix && y == iy)
662
0
            return 1;
663
0
    }
664
0
    return 0;
665
0
}
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
0
{
679
0
l_int32  i, j, n1, n2, x1, y1, x2, y2;
680
681
0
    if (!pta1)
682
0
        return ERROR_INT("pta1 not defined", __func__, 0);
683
0
    if (!pta2)
684
0
        return ERROR_INT("pta2 not defined", __func__, 0);
685
686
0
    n1 = ptaGetCount(pta1);
687
0
    n2 = ptaGetCount(pta2);
688
0
    for (i = 0; i < n1; i++) {
689
0
        ptaGetIPt(pta1, i, &x1, &y1);
690
0
        for (j = 0; j < n2; j++) {
691
0
            ptaGetIPt(pta2, i, &x2, &y2);
692
0
            if (x1 == x2 && y1 == y2)
693
0
                return 1;
694
0
        }
695
0
    }
696
697
0
    return 0;
698
0
}
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
0
{
721
0
l_int32  n, i, x, y;
722
0
PTA     *ptad;
723
724
0
    if (!ptas)
725
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
726
0
    n = ptaGetCount(ptas);
727
0
    ptad = ptaCreate(n);
728
0
    for (i = 0; i < n; i++) {
729
0
        ptaGetIPt(ptas, i, &x, &y);
730
0
        x = (l_int32)(scalex * (x + shiftx) + 0.5);
731
0
        y = (l_int32)(scaley * (y + shifty) + 0.5);
732
0
        ptaAddPt(ptad, x, y);
733
0
    }
734
735
0
    return ptad;
736
0
}
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
0
{
758
0
l_int32    i, n;
759
0
l_float32  sum, x1, y1, x2, y2, xp1, yp1, xp2, yp2;
760
761
0
    if (!pinside)
762
0
        return ERROR_INT("&inside not defined", __func__, 1);
763
0
    *pinside = 0;
764
0
    if (!pta)
765
0
        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
0
    n = ptaGetCount(pta);
770
0
    sum = 0.0;
771
0
    for (i = 0; i < n; i++) {
772
0
        ptaGetPt(pta, i, &xp1, &yp1);
773
0
        ptaGetPt(pta, (i + 1) % n, &xp2, &yp2);
774
0
        x1 = xp1 - x;
775
0
        y1 = yp1 - y;
776
0
        x2 = xp2 - x;
777
0
        y2 = yp2 - y;
778
0
        sum += l_angleBetweenVectors(x1, y1, x2, y2);
779
0
    }
780
781
0
    if (L_ABS(sum) > M_PI)
782
0
        *pinside = 1;
783
0
    return 0;
784
0
}
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
0
{
808
0
l_float64  ang;
809
810
0
    ang = atan2(y2, x2) - atan2(y1, x1);
811
0
    if (ang > M_PI) ang -= 2.0 * M_PI;
812
0
    if (ang < -M_PI) ang += 2.0 * M_PI;
813
0
    return ang;
814
0
}
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
0
{
841
0
l_int32    i, n;
842
0
l_float32  x0, y0, x1, y1, x2, y2;
843
0
l_float64  cprod;
844
845
0
    if (!pisconvex)
846
0
        return ERROR_INT("&isconvex not defined", __func__, 1);
847
0
    *pisconvex = 0;
848
0
    if (!pta)
849
0
        return ERROR_INT("pta not defined", __func__, 1);
850
0
    if ((n = ptaGetCount(pta)) < 3)
851
0
        return ERROR_INT("pta has < 3 pts", __func__, 1);
852
853
0
    for (i = 0; i < n; i++) {
854
0
        ptaGetPt(pta, i, &x0, &y0);
855
0
        ptaGetPt(pta, (i + 1) % n, &x1, &y1);
856
0
        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
0
        cprod = (x2 - x0) * (y1 - y0) - (x1 - x0) * (y2 - y0);
864
0
        if (cprod < -0.0001)  /* small delta for float accuracy; test fails */
865
0
            return 0;
866
0
    }
867
0
    *pisconvex = 1;
868
0
    return 0;
869
0
}
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
0
{
893
0
l_int32    i, n;
894
0
l_float32  x, y, xmin, ymin, xmax, ymax;
895
896
0
    if (pxmin) *pxmin = -1.0;
897
0
    if (pymin) *pymin = -1.0;
898
0
    if (pxmax) *pxmax = -1.0;
899
0
    if (pymax) *pymax = -1.0;
900
0
    if (!pta)
901
0
        return ERROR_INT("pta not defined", __func__, 1);
902
0
    if (!pxmin && !pxmax && !pymin && !pymax)
903
0
        return ERROR_INT("no output requested", __func__, 1);
904
0
    if ((n = ptaGetCount(pta)) == 0) {
905
0
        L_WARNING("pta is empty\n", __func__);
906
0
        return 0;
907
0
    }
908
909
0
    xmin = ymin = 1.0e20f;
910
0
    xmax = ymax = -1.0e20f;
911
0
    for (i = 0; i < n; i++) {
912
0
        ptaGetPt(pta, i, &x, &y);
913
0
        if (x < xmin) xmin = x;
914
0
        if (y < ymin) ymin = y;
915
0
        if (x > xmax) xmax = x;
916
0
        if (y > ymax) ymax = y;
917
0
    }
918
0
    if (pxmin) *pxmin = xmin;
919
0
    if (pymin) *pymin = ymin;
920
0
    if (pxmax) *pxmax = xmax;
921
0
    if (pymax) *pymax = ymax;
922
0
    return 0;
923
0
}
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
0
{
944
0
l_int32    i, n;
945
0
l_float32  x, y;
946
0
PTA       *ptad;
947
948
0
    if (!ptas)
949
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
950
0
    if (ptaGetCount(ptas) == 0) {
951
0
        L_WARNING("ptas is empty\n", __func__);
952
0
        return ptaCopy(ptas);
953
0
    }
954
0
    if (type != L_SELECT_XVAL && type != L_SELECT_YVAL &&
955
0
        type != L_SELECT_IF_EITHER && type != L_SELECT_IF_BOTH)
956
0
        return (PTA *)ERROR_PTR("invalid type", __func__, NULL);
957
0
    if (relation != L_SELECT_IF_LT && relation != L_SELECT_IF_GT &&
958
0
        relation != L_SELECT_IF_LTE && relation != L_SELECT_IF_GTE)
959
0
        return (PTA *)ERROR_PTR("invalid relation", __func__, NULL);
960
961
0
    n = ptaGetCount(ptas);
962
0
    ptad = ptaCreate(n);
963
0
    for (i = 0; i < n; i++) {
964
0
        ptaGetPt(ptas, i, &x, &y);
965
0
        if (type == L_SELECT_XVAL) {
966
0
            if ((relation == L_SELECT_IF_LT && x < xth) ||
967
0
                (relation == L_SELECT_IF_GT && x > xth) ||
968
0
                (relation == L_SELECT_IF_LTE && x <= xth) ||
969
0
                (relation == L_SELECT_IF_GTE && x >= xth))
970
0
                ptaAddPt(ptad, x, y);
971
0
        } else if (type == L_SELECT_YVAL) {
972
0
            if ((relation == L_SELECT_IF_LT && y < yth) ||
973
0
                (relation == L_SELECT_IF_GT && y > yth) ||
974
0
                (relation == L_SELECT_IF_LTE && y <= yth) ||
975
0
                (relation == L_SELECT_IF_GTE && y >= yth))
976
0
                ptaAddPt(ptad, x, y);
977
0
        } else if (type == L_SELECT_IF_EITHER) {
978
0
            if (((relation == L_SELECT_IF_LT) && (x < xth || y < yth)) ||
979
0
                ((relation == L_SELECT_IF_GT) && (x > xth || y > yth)) ||
980
0
                ((relation == L_SELECT_IF_LTE) && (x <= xth || y <= yth)) ||
981
0
                ((relation == L_SELECT_IF_GTE) && (x >= xth || y >= yth)))
982
0
                ptaAddPt(ptad, x, y);
983
0
        } else {  /* L_SELECT_IF_BOTH */
984
0
            if (((relation == L_SELECT_IF_LT) && (x < xth && y < yth)) ||
985
0
                ((relation == L_SELECT_IF_GT) && (x > xth && y > yth)) ||
986
0
                ((relation == L_SELECT_IF_LTE) && (x <= xth && y <= yth)) ||
987
0
                ((relation == L_SELECT_IF_GTE) && (x >= xth && y >= yth)))
988
0
                ptaAddPt(ptad, x, y);
989
0
        }
990
0
    }
991
992
0
    return ptad;
993
0
}
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
0
{
1007
0
l_int32   i, n, x, y;
1008
0
l_uint32  val;
1009
0
PTA      *ptad;
1010
1011
0
    if (!ptas)
1012
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
1013
0
    if (!pixm || pixGetDepth(pixm) != 1)
1014
0
        return (PTA *)ERROR_PTR("pixm undefined or not 1 bpp", __func__, NULL);
1015
0
    if (ptaGetCount(ptas) == 0) {
1016
0
        L_INFO("ptas is empty\n", __func__);
1017
0
        return ptaCopy(ptas);
1018
0
    }
1019
1020
0
    n = ptaGetCount(ptas);
1021
0
    ptad = ptaCreate(n);
1022
0
    for (i = 0; i < n; i++) {
1023
0
        ptaGetIPt(ptas, i, &x, &y);
1024
0
        pixGetPixel(pixm, x, y, &val);
1025
0
        if (val == 1)
1026
0
            ptaAddPt(ptad, x, y);
1027
0
    }
1028
0
    return ptad;
1029
0
}
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
0
{
1074
0
l_int32     n, i;
1075
0
l_float32   a, b, factor, sx, sy, sxx, sxy, val;
1076
0
l_float32  *xa, *ya;
1077
1078
0
    if (pa) *pa = 0.0;
1079
0
    if (pb) *pb = 0.0;
1080
0
    if (pnafit) *pnafit = NULL;
1081
0
    if (!pa && !pb && !pnafit)
1082
0
        return ERROR_INT("no output requested", __func__, 1);
1083
0
    if (!pta)
1084
0
        return ERROR_INT("pta not defined", __func__, 1);
1085
0
    if ((n = ptaGetCount(pta)) < 2)
1086
0
        return ERROR_INT("less than 2 pts found", __func__, 1);
1087
1088
0
    xa = pta->x;  /* not a copy */
1089
0
    ya = pta->y;  /* not a copy */
1090
0
    sx = sy = sxx = sxy = 0.;
1091
0
    if (pa && pb) {  /* general line */
1092
0
        for (i = 0; i < n; i++) {
1093
0
            sx += xa[i];
1094
0
            sy += ya[i];
1095
0
            sxx += xa[i] * xa[i];
1096
0
            sxy += xa[i] * ya[i];
1097
0
        }
1098
0
        factor = n * sxx - sx * sx;
1099
0
        if (factor == 0.0)
1100
0
            return ERROR_INT("no solution found", __func__, 1);
1101
0
        factor = 1.f / factor;
1102
1103
0
        a = factor * ((l_float32)n * sxy - sx * sy);
1104
0
        b = factor * (sxx * sy - sx * sxy);
1105
0
    } else if (pa) {  /* b = 0; line through origin */
1106
0
        for (i = 0; i < n; i++) {
1107
0
            sxx += xa[i] * xa[i];
1108
0
            sxy += xa[i] * ya[i];
1109
0
        }
1110
0
        if (sxx == 0.0)
1111
0
            return ERROR_INT("no solution found", __func__, 1);
1112
0
        a = sxy / sxx;
1113
0
        b = 0.0;
1114
0
    } else {  /* a = 0; horizontal line */
1115
0
        for (i = 0; i < n; i++)
1116
0
            sy += ya[i];
1117
0
        a = 0.0;
1118
0
        b = sy / (l_float32)n;
1119
0
    }
1120
1121
0
    if (pnafit) {
1122
0
        *pnafit = numaCreate(n);
1123
0
        for (i = 0; i < n; i++) {
1124
0
            val = a * xa[i] + b;
1125
0
            numaAddNumber(*pnafit, val);
1126
0
        }
1127
0
    }
1128
1129
0
    if (pa) *pa = a;
1130
0
    if (pb) *pb = b;
1131
0
    return 0;
1132
0
}
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
0
{
1174
0
l_int32     n, i, ret;
1175
0
l_float32   x, y, sx, sy, sx2, sx3, sx4, sxy, sx2y;
1176
0
l_float32  *xa, *ya;
1177
0
l_float32  *f[3];
1178
0
l_float32   g[3];
1179
1180
0
    if (pa) *pa = 0.0;
1181
0
    if (pb) *pb = 0.0;
1182
0
    if (pc) *pc = 0.0;
1183
0
    if (pnafit) *pnafit = NULL;
1184
0
    if (!pa && !pb && !pc && !pnafit)
1185
0
        return ERROR_INT("no output requested", __func__, 1);
1186
0
    if (!pta)
1187
0
        return ERROR_INT("pta not defined", __func__, 1);
1188
0
    if ((n = ptaGetCount(pta)) < 3)
1189
0
        return ERROR_INT("less than 3 pts found", __func__, 1);
1190
1191
0
    xa = pta->x;  /* not a copy */
1192
0
    ya = pta->y;  /* not a copy */
1193
0
    sx = sy = sx2 = sx3 = sx4 = sxy = sx2y = 0.;
1194
0
    for (i = 0; i < n; i++) {
1195
0
        x = xa[i];
1196
0
        y = ya[i];
1197
0
        sx += x;
1198
0
        sy += y;
1199
0
        sx2 += x * x;
1200
0
        sx3 += x * x * x;
1201
0
        sx4 += x * x * x * x;
1202
0
        sxy += x * y;
1203
0
        sx2y += x * x * y;
1204
0
    }
1205
1206
0
    for (i = 0; i < 3; i++)
1207
0
        f[i] = (l_float32 *)LEPT_CALLOC(3, sizeof(l_float32));
1208
0
    f[0][0] = sx4;
1209
0
    f[0][1] = sx3;
1210
0
    f[0][2] = sx2;
1211
0
    f[1][0] = sx3;
1212
0
    f[1][1] = sx2;
1213
0
    f[1][2] = sx;
1214
0
    f[2][0] = sx2;
1215
0
    f[2][1] = sx;
1216
0
    f[2][2] = n;
1217
0
    g[0] = sx2y;
1218
0
    g[1] = sxy;
1219
0
    g[2] = sy;
1220
1221
        /* Solve for the unknowns, also putting f-inverse into f */
1222
0
    ret = gaussjordan(f, g, 3);
1223
0
    for (i = 0; i < 3; i++)
1224
0
        LEPT_FREE(f[i]);
1225
0
    if (ret)
1226
0
        return ERROR_INT("quadratic solution failed", __func__, 1);
1227
1228
0
    if (pa) *pa = g[0];
1229
0
    if (pb) *pb = g[1];
1230
0
    if (pc) *pc = g[2];
1231
0
    if (pnafit) {
1232
0
        *pnafit = numaCreate(n);
1233
0
        for (i = 0; i < n; i++) {
1234
0
            x = xa[i];
1235
0
            y = g[0] * x * x + g[1] * x + g[2];
1236
0
            numaAddNumber(*pnafit, y);
1237
0
        }
1238
0
    }
1239
0
    return 0;
1240
0
}
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
0
{
1286
0
l_int32     n, i, ret;
1287
0
l_float32   x, y, sx, sy, sx2, sx3, sx4, sx5, sx6, sxy, sx2y, sx3y;
1288
0
l_float32  *xa, *ya;
1289
0
l_float32  *f[4];
1290
0
l_float32   g[4];
1291
1292
0
    if (pa) *pa = 0.0;
1293
0
    if (pb) *pb = 0.0;
1294
0
    if (pc) *pc = 0.0;
1295
0
    if (pd) *pd = 0.0;
1296
0
    if (pnafit) *pnafit = NULL;
1297
0
    if (!pa && !pb && !pc && !pd && !pnafit)
1298
0
        return ERROR_INT("no output requested", __func__, 1);
1299
0
    if (!pta)
1300
0
        return ERROR_INT("pta not defined", __func__, 1);
1301
0
    if ((n = ptaGetCount(pta)) < 4)
1302
0
        return ERROR_INT("less than 4 pts found", __func__, 1);
1303
1304
0
    xa = pta->x;  /* not a copy */
1305
0
    ya = pta->y;  /* not a copy */
1306
0
    sx = sy = sx2 = sx3 = sx4 = sx5 = sx6 = sxy = sx2y = sx3y = 0.;
1307
0
    for (i = 0; i < n; i++) {
1308
0
        x = xa[i];
1309
0
        y = ya[i];
1310
0
        sx += x;
1311
0
        sy += y;
1312
0
        sx2 += x * x;
1313
0
        sx3 += x * x * x;
1314
0
        sx4 += x * x * x * x;
1315
0
        sx5 += x * x * x * x * x;
1316
0
        sx6 += x * x * x * x * x * x;
1317
0
        sxy += x * y;
1318
0
        sx2y += x * x * y;
1319
0
        sx3y += x * x * x * y;
1320
0
    }
1321
1322
0
    for (i = 0; i < 4; i++)
1323
0
        f[i] = (l_float32 *)LEPT_CALLOC(4, sizeof(l_float32));
1324
0
    f[0][0] = sx6;
1325
0
    f[0][1] = sx5;
1326
0
    f[0][2] = sx4;
1327
0
    f[0][3] = sx3;
1328
0
    f[1][0] = sx5;
1329
0
    f[1][1] = sx4;
1330
0
    f[1][2] = sx3;
1331
0
    f[1][3] = sx2;
1332
0
    f[2][0] = sx4;
1333
0
    f[2][1] = sx3;
1334
0
    f[2][2] = sx2;
1335
0
    f[2][3] = sx;
1336
0
    f[3][0] = sx3;
1337
0
    f[3][1] = sx2;
1338
0
    f[3][2] = sx;
1339
0
    f[3][3] = n;
1340
0
    g[0] = sx3y;
1341
0
    g[1] = sx2y;
1342
0
    g[2] = sxy;
1343
0
    g[3] = sy;
1344
1345
        /* Solve for the unknowns, also putting f-inverse into f */
1346
0
    ret = gaussjordan(f, g, 4);
1347
0
    for (i = 0; i < 4; i++)
1348
0
        LEPT_FREE(f[i]);
1349
0
    if (ret)
1350
0
        return ERROR_INT("cubic solution failed", __func__, 1);
1351
1352
0
    if (pa) *pa = g[0];
1353
0
    if (pb) *pb = g[1];
1354
0
    if (pc) *pc = g[2];
1355
0
    if (pd) *pd = g[3];
1356
0
    if (pnafit) {
1357
0
        *pnafit = numaCreate(n);
1358
0
        for (i = 0; i < n; i++) {
1359
0
            x = xa[i];
1360
0
            y = g[0] * x * x * x + g[1] * x * x + g[2] * x + g[3];
1361
0
            numaAddNumber(*pnafit, y);
1362
0
        }
1363
0
    }
1364
0
    return 0;
1365
0
}
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
0
{
1415
0
l_int32     n, i, ret;
1416
0
l_float32   x, y, sx, sy, sx2, sx3, sx4, sx5, sx6, sx7, sx8;
1417
0
l_float32   sxy, sx2y, sx3y, sx4y;
1418
0
l_float32  *xa, *ya;
1419
0
l_float32  *f[5];
1420
0
l_float32   g[5];
1421
1422
0
    if (pa) *pa = 0.0;
1423
0
    if (pb) *pb = 0.0;
1424
0
    if (pc) *pc = 0.0;
1425
0
    if (pd) *pd = 0.0;
1426
0
    if (pe) *pe = 0.0;
1427
0
    if (pnafit) *pnafit = NULL;
1428
0
    if (!pa && !pb && !pc && !pd && !pe && !pnafit)
1429
0
        return ERROR_INT("no output requested", __func__, 1);
1430
0
    if (!pta)
1431
0
        return ERROR_INT("pta not defined", __func__, 1);
1432
0
    if ((n = ptaGetCount(pta)) < 5)
1433
0
        return ERROR_INT("less than 5 pts found", __func__, 1);
1434
1435
0
    xa = pta->x;  /* not a copy */
1436
0
    ya = pta->y;  /* not a copy */
1437
0
    sx = sy = sx2 = sx3 = sx4 = sx5 = sx6 = sx7 = sx8 = 0;
1438
0
    sxy = sx2y = sx3y = sx4y = 0.;
1439
0
    for (i = 0; i < n; i++) {
1440
0
        x = xa[i];
1441
0
        y = ya[i];
1442
0
        sx += x;
1443
0
        sy += y;
1444
0
        sx2 += x * x;
1445
0
        sx3 += x * x * x;
1446
0
        sx4 += x * x * x * x;
1447
0
        sx5 += x * x * x * x * x;
1448
0
        sx6 += x * x * x * x * x * x;
1449
0
        sx7 += x * x * x * x * x * x * x;
1450
0
        sx8 += x * x * x * x * x * x * x * x;
1451
0
        sxy += x * y;
1452
0
        sx2y += x * x * y;
1453
0
        sx3y += x * x * x * y;
1454
0
        sx4y += x * x * x * x * y;
1455
0
    }
1456
1457
0
    for (i = 0; i < 5; i++)
1458
0
        f[i] = (l_float32 *)LEPT_CALLOC(5, sizeof(l_float32));
1459
0
    f[0][0] = sx8;
1460
0
    f[0][1] = sx7;
1461
0
    f[0][2] = sx6;
1462
0
    f[0][3] = sx5;
1463
0
    f[0][4] = sx4;
1464
0
    f[1][0] = sx7;
1465
0
    f[1][1] = sx6;
1466
0
    f[1][2] = sx5;
1467
0
    f[1][3] = sx4;
1468
0
    f[1][4] = sx3;
1469
0
    f[2][0] = sx6;
1470
0
    f[2][1] = sx5;
1471
0
    f[2][2] = sx4;
1472
0
    f[2][3] = sx3;
1473
0
    f[2][4] = sx2;
1474
0
    f[3][0] = sx5;
1475
0
    f[3][1] = sx4;
1476
0
    f[3][2] = sx3;
1477
0
    f[3][3] = sx2;
1478
0
    f[3][4] = sx;
1479
0
    f[4][0] = sx4;
1480
0
    f[4][1] = sx3;
1481
0
    f[4][2] = sx2;
1482
0
    f[4][3] = sx;
1483
0
    f[4][4] = n;
1484
0
    g[0] = sx4y;
1485
0
    g[1] = sx3y;
1486
0
    g[2] = sx2y;
1487
0
    g[3] = sxy;
1488
0
    g[4] = sy;
1489
1490
        /* Solve for the unknowns, also putting f-inverse into f */
1491
0
    ret = gaussjordan(f, g, 5);
1492
0
    for (i = 0; i < 5; i++)
1493
0
        LEPT_FREE(f[i]);
1494
0
    if (ret)
1495
0
        return ERROR_INT("quartic solution failed", __func__, 1);
1496
1497
0
    if (pa) *pa = g[0];
1498
0
    if (pb) *pb = g[1];
1499
0
    if (pc) *pc = g[2];
1500
0
    if (pd) *pd = g[3];
1501
0
    if (pe) *pe = g[4];
1502
0
    if (pnafit) {
1503
0
        *pnafit = numaCreate(n);
1504
0
        for (i = 0; i < n; i++) {
1505
0
            x = xa[i];
1506
0
            y = g[0] * x * x * x * x + g[1] * x * x * x + g[2] * x * x
1507
0
                 + g[3] * x + g[4];
1508
0
            numaAddNumber(*pnafit, y);
1509
0
        }
1510
0
    }
1511
0
    return 0;
1512
0
}
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
0
{
1549
0
l_int32    n, i, ret;
1550
0
l_float32  x, y, yf, val, mederr;
1551
0
NUMA      *nafit, *naerror;
1552
0
PTA       *ptad;
1553
1554
0
    if (pptad) *pptad = NULL;
1555
0
    if (pa) *pa = 0.0;
1556
0
    if (pb) *pb = 0.0;
1557
0
    if (pmederr) *pmederr = 0.0;
1558
0
    if (pnafit) *pnafit = NULL;
1559
0
    if (!pptad && !pa && !pb && !pnafit)
1560
0
        return ERROR_INT("no output requested", __func__, 1);
1561
0
    if (!pta)
1562
0
        return ERROR_INT("pta not defined", __func__, 1);
1563
0
    if (factor <= 0.0)
1564
0
        return ERROR_INT("factor must be > 0.0", __func__, 1);
1565
0
    if ((n = ptaGetCount(pta)) < 3)
1566
0
        return ERROR_INT("less than 2 pts found", __func__, 1);
1567
1568
0
    if (ptaGetLinearLSF(pta, pa, pb, &nafit) != 0)
1569
0
        return ERROR_INT("error in linear LSF", __func__, 1);
1570
1571
        /* Get the median error */
1572
0
    naerror = numaCreate(n);
1573
0
    for (i = 0; i < n; i++) {
1574
0
        ptaGetPt(pta, i, &x, &y);
1575
0
        numaGetFValue(nafit, i, &yf);
1576
0
        numaAddNumber(naerror, L_ABS(y - yf));
1577
0
    }
1578
0
    numaGetMedian(naerror, &mederr);
1579
0
    if (pmederr) *pmederr = mederr;
1580
0
    numaDestroy(&nafit);
1581
1582
        /* Remove outliers */
1583
0
    ptad = ptaCreate(n);
1584
0
    for (i = 0; i < n; i++) {
1585
0
        ptaGetPt(pta, i, &x, &y);
1586
0
        numaGetFValue(naerror, i, &val);
1587
0
        if (val <= factor * mederr)  /* <= in case mederr = 0 */
1588
0
            ptaAddPt(ptad, x, y);
1589
0
    }
1590
0
    numaDestroy(&naerror);
1591
1592
       /* Do LSF again */
1593
0
    ret = ptaGetLinearLSF(ptad, pa, pb, pnafit);
1594
0
    if (pptad)
1595
0
        *pptad = ptad;
1596
0
    else
1597
0
        ptaDestroy(&ptad);
1598
1599
0
    return ret;
1600
0
}
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
0
{
1635
0
l_int32    n, i, ret;
1636
0
l_float32  x, y, yf, val, mederr;
1637
0
NUMA      *nafit, *naerror;
1638
0
PTA       *ptad;
1639
1640
0
    if (pptad) *pptad = NULL;
1641
0
    if (pa) *pa = 0.0;
1642
0
    if (pb) *pb = 0.0;
1643
0
    if (pc) *pc = 0.0;
1644
0
    if (pmederr) *pmederr = 0.0;
1645
0
    if (pnafit) *pnafit = NULL;
1646
0
    if (!pptad && !pa && !pb && !pc && !pnafit)
1647
0
        return ERROR_INT("no output requested", __func__, 1);
1648
0
    if (factor <= 0.0)
1649
0
        return ERROR_INT("factor must be > 0.0", __func__, 1);
1650
0
    if (!pta)
1651
0
        return ERROR_INT("pta not defined", __func__, 1);
1652
0
    if ((n = ptaGetCount(pta)) < 3)
1653
0
        return ERROR_INT("less than 3 pts found", __func__, 1);
1654
1655
0
    if (ptaGetQuadraticLSF(pta, NULL, NULL, NULL, &nafit) != 0)
1656
0
        return ERROR_INT("error in quadratic LSF", __func__, 1);
1657
1658
        /* Get the median error */
1659
0
    naerror = numaCreate(n);
1660
0
    for (i = 0; i < n; i++) {
1661
0
        ptaGetPt(pta, i, &x, &y);
1662
0
        numaGetFValue(nafit, i, &yf);
1663
0
        numaAddNumber(naerror, L_ABS(y - yf));
1664
0
    }
1665
0
    numaGetMedian(naerror, &mederr);
1666
0
    if (pmederr) *pmederr = mederr;
1667
0
    numaDestroy(&nafit);
1668
1669
        /* Remove outliers */
1670
0
    ptad = ptaCreate(n);
1671
0
    for (i = 0; i < n; i++) {
1672
0
        ptaGetPt(pta, i, &x, &y);
1673
0
        numaGetFValue(naerror, i, &val);
1674
0
        if (val <= factor * mederr)  /* <= in case mederr = 0 */
1675
0
            ptaAddPt(ptad, x, y);
1676
0
    }
1677
0
    numaDestroy(&naerror);
1678
0
    n = ptaGetCount(ptad);
1679
0
    if ((n = ptaGetCount(ptad)) < 3) {
1680
0
        ptaDestroy(&ptad);
1681
0
        return ERROR_INT("less than 3 pts found", __func__, 1);
1682
0
    }
1683
1684
       /* Do LSF again */
1685
0
    ret = ptaGetQuadraticLSF(ptad, pa, pb, pc, pnafit);
1686
0
    if (pptad)
1687
0
        *pptad = ptad;
1688
0
    else
1689
0
        ptaDestroy(&ptad);
1690
1691
0
    return ret;
1692
0
}
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
0
{
1709
0
    if (!py)
1710
0
        return ERROR_INT("&y not defined", __func__, 1);
1711
1712
0
    *py = a * x + b;
1713
0
    return 0;
1714
0
}
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
0
{
1732
0
    if (!py)
1733
0
        return ERROR_INT("&y not defined", __func__, 1);
1734
1735
0
    *py = a * x * x + b * x + c;
1736
0
    return 0;
1737
0
}
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
0
{
1756
0
    if (!py)
1757
0
        return ERROR_INT("&y not defined", __func__, 1);
1758
1759
0
    *py = a * x * x * x + b * x * x + c * x + d;
1760
0
    return 0;
1761
0
}
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
0
{
1781
0
l_float32  x2;
1782
1783
0
    if (!py)
1784
0
        return ERROR_INT("&y not defined", __func__, 1);
1785
1786
0
    x2 = x * x;
1787
0
    *py = a * x2 * x2 + b * x2 * x + c * x2 + d * x + e;
1788
0
    return 0;
1789
0
}
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
0
{
1817
0
char            buffer[128];
1818
0
char           *rtitle, *gtitle, *btitle;
1819
0
static l_int32  count = 0;  /* require separate temp files for each call */
1820
0
l_int32         i, x, y, d, w, h, npts, rval, gval, bval;
1821
0
l_uint32        val;
1822
0
NUMA           *na, *nar, *nag, *nab;
1823
0
PIX            *pixt;
1824
1825
0
    lept_mkdir("lept/plot");
1826
1827
0
    if (!pixs)
1828
0
        return ERROR_INT("pixs not defined", __func__, 1);
1829
0
    if (!pta)
1830
0
        return ERROR_INT("pta not defined", __func__, 1);
1831
0
    if (outformat != GPLOT_PNG && outformat != GPLOT_PS &&
1832
0
        outformat != GPLOT_EPS && outformat != GPLOT_LATEX) {
1833
0
        L_WARNING("outformat invalid; using GPLOT_PNG\n", __func__);
1834
0
        outformat = GPLOT_PNG;
1835
0
    }
1836
1837
0
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
1838
0
    d = pixGetDepth(pixt);
1839
0
    w = pixGetWidth(pixt);
1840
0
    h = pixGetHeight(pixt);
1841
0
    npts = ptaGetCount(pta);
1842
0
    if (d == 32) {
1843
0
        nar = numaCreate(npts);
1844
0
        nag = numaCreate(npts);
1845
0
        nab = numaCreate(npts);
1846
0
        for (i = 0; i < npts; i++) {
1847
0
            ptaGetIPt(pta, i, &x, &y);
1848
0
            if (x < 0 || x >= w)
1849
0
                continue;
1850
0
            if (y < 0 || y >= h)
1851
0
                continue;
1852
0
            pixGetPixel(pixt, x, y, &val);
1853
0
            rval = GET_DATA_BYTE(&val, COLOR_RED);
1854
0
            gval = GET_DATA_BYTE(&val, COLOR_GREEN);
1855
0
            bval = GET_DATA_BYTE(&val, COLOR_BLUE);
1856
0
            numaAddNumber(nar, rval);
1857
0
            numaAddNumber(nag, gval);
1858
0
            numaAddNumber(nab, bval);
1859
0
        }
1860
1861
0
        snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
1862
0
        rtitle = stringJoin("Red: ", title);
1863
0
        gplotSimple1(nar, outformat, buffer, rtitle);
1864
0
        snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
1865
0
        gtitle = stringJoin("Green: ", title);
1866
0
        gplotSimple1(nag, outformat, buffer, gtitle);
1867
0
        snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
1868
0
        btitle = stringJoin("Blue: ", title);
1869
0
        gplotSimple1(nab, outformat, buffer, btitle);
1870
0
        numaDestroy(&nar);
1871
0
        numaDestroy(&nag);
1872
0
        numaDestroy(&nab);
1873
0
        LEPT_FREE(rtitle);
1874
0
        LEPT_FREE(gtitle);
1875
0
        LEPT_FREE(btitle);
1876
0
    } else {
1877
0
        na = numaCreate(npts);
1878
0
        for (i = 0; i < npts; i++) {
1879
0
            ptaGetIPt(pta, i, &x, &y);
1880
0
            if (x < 0 || x >= w)
1881
0
                continue;
1882
0
            if (y < 0 || y >= h)
1883
0
                continue;
1884
0
            pixGetPixel(pixt, x, y, &val);
1885
0
            numaAddNumber(na, (l_float32)val);
1886
0
        }
1887
1888
0
        snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
1889
0
        gplotSimple1(na, outformat, buffer, title);
1890
0
        numaDestroy(&na);
1891
0
    }
1892
0
    pixDestroy(&pixt);
1893
0
    return 0;
1894
0
}
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
0
{
1914
0
l_int32    i, j, w, h, wpl, xstart, xend, ystart, yend, bw, bh;
1915
0
l_uint32  *data, *line;
1916
0
PTA       *pta;
1917
1918
0
    if (!pixs || (pixGetDepth(pixs) != 1))
1919
0
        return (PTA *)ERROR_PTR("pixs undefined or not 1 bpp", __func__, NULL);
1920
1921
0
    pixGetDimensions(pixs, &w, &h, NULL);
1922
0
    data = pixGetData(pixs);
1923
0
    wpl = pixGetWpl(pixs);
1924
0
    xstart = ystart = 0;
1925
0
    xend = w - 1;
1926
0
    yend = h - 1;
1927
0
    if (box) {
1928
0
        boxGetGeometry(box, &xstart, &ystart, &bw, &bh);
1929
0
        xend = xstart + bw - 1;
1930
0
        yend = ystart + bh - 1;
1931
0
    }
1932
1933
0
    if ((pta = ptaCreate(0)) == NULL)
1934
0
        return (PTA *)ERROR_PTR("pta not made", __func__, NULL);
1935
0
    for (i = ystart; i <= yend; i++) {
1936
0
        line = data + i * wpl;
1937
0
        for (j = xstart; j <= xend; j++) {
1938
0
            if (GET_DATA_BIT(line, j))
1939
0
                ptaAddPt(pta, j, i);
1940
0
        }
1941
0
    }
1942
1943
0
    return pta;
1944
0
}
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
0
{
1966
0
l_int32  n, i, x, y;
1967
0
PIX     *pix;
1968
1969
0
    if (!pta)
1970
0
        return (PIX *)ERROR_PTR("pta not defined", __func__, NULL);
1971
1972
0
    if ((pix = pixCreate(w, h, 1)) == NULL)
1973
0
        return (PIX *)ERROR_PTR("pix not made", __func__, NULL);
1974
0
    n = ptaGetCount(pta);
1975
0
    for (i = 0; i < n; i++) {
1976
0
        ptaGetIPt(pta, i, &x, &y);
1977
0
        if (x < 0 || x >= w || y < 0 || y >= h)
1978
0
            continue;
1979
0
        pixSetPixel(pix, x, y, 1);
1980
0
    }
1981
1982
0
    return pix;
1983
0
}
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
0
{
2004
0
PIX  *pixt;
2005
0
PTA  *pta;
2006
2007
0
    if (!pixs || (pixGetDepth(pixs) != 1))
2008
0
        return (PTA *)ERROR_PTR("pixs undefined or not 1 bpp", __func__, NULL);
2009
0
    if (type != L_BOUNDARY_FG && type != L_BOUNDARY_BG)
2010
0
        return (PTA *)ERROR_PTR("invalid type", __func__, NULL);
2011
2012
0
    if (type == L_BOUNDARY_FG)
2013
0
        pixt = pixMorphSequence(pixs, "e3.3", 0);
2014
0
    else
2015
0
        pixt = pixMorphSequence(pixs, "d3.3", 0);
2016
0
    pixXor(pixt, pixt, pixs);
2017
0
    pta = ptaGetPixelsFromPix(pixt, NULL);
2018
2019
0
    pixDestroy(&pixt);
2020
0
    return pta;
2021
0
}
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
0
{
2054
0
l_int32  i, n, w, h, x, y, bw, bh, left, right, top, bot;
2055
0
BOXA    *boxa;
2056
0
PIX     *pixt1, *pixt2;
2057
0
PIXA    *pixa;
2058
0
PTA     *pta1, *pta2;
2059
0
PTAA    *ptaa;
2060
2061
0
    if (pboxa) *pboxa = NULL;
2062
0
    if (ppixa) *ppixa = NULL;
2063
0
    if (!pixs || (pixGetDepth(pixs) != 1))
2064
0
        return (PTAA *)ERROR_PTR("pixs undefined or not 1 bpp", __func__, NULL);
2065
0
    if (type != L_BOUNDARY_FG && type != L_BOUNDARY_BG)
2066
0
        return (PTAA *)ERROR_PTR("invalid type", __func__, NULL);
2067
0
    if (connectivity != 4 && connectivity != 8)
2068
0
        return (PTAA *)ERROR_PTR("connectivity not 4 or 8", __func__, NULL);
2069
2070
0
    pixGetDimensions(pixs, &w, &h, NULL);
2071
0
    boxa = pixConnComp(pixs, &pixa, connectivity);
2072
0
    n = boxaGetCount(boxa);
2073
0
    ptaa = ptaaCreate(0);
2074
0
    for (i = 0; i < n; i++) {
2075
0
        pixt1 = pixaGetPix(pixa, i, L_CLONE);
2076
0
        boxaGetBoxGeometry(boxa, i, &x, &y, &bw, &bh);
2077
0
        left = right = top = bot = 0;
2078
0
        if (type == L_BOUNDARY_BG) {
2079
0
            if (x > 0) left = 1;
2080
0
            if (y > 0) top = 1;
2081
0
            if (x + bw < w) right = 1;
2082
0
            if (y + bh < h) bot = 1;
2083
0
            pixt2 = pixAddBorderGeneral(pixt1, left, right, top, bot, 0);
2084
0
        } else {
2085
0
            pixt2 = pixClone(pixt1);
2086
0
        }
2087
0
        pta1 = ptaGetBoundaryPixels(pixt2, type);
2088
0
        pta2 = ptaTransform(pta1, x - left, y - top, 1.0, 1.0);
2089
0
        ptaaAddPta(ptaa, pta2, L_INSERT);
2090
0
        ptaDestroy(&pta1);
2091
0
        pixDestroy(&pixt1);
2092
0
        pixDestroy(&pixt2);
2093
0
    }
2094
2095
0
    if (pboxa)
2096
0
        *pboxa = boxa;
2097
0
    else
2098
0
        boxaDestroy(&boxa);
2099
0
    if (ppixa)
2100
0
        *ppixa = pixa;
2101
0
    else
2102
0
        pixaDestroy(&pixa);
2103
0
    return ptaa;
2104
0
}
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
0
{
2132
0
l_int32    wpl, index, i, j, w, h;
2133
0
l_uint32   maxval;
2134
0
l_uint32  *data, *line;
2135
0
PTA       *pta;
2136
0
PTAA      *ptaa;
2137
2138
0
    if (pncc) *pncc = 0;
2139
0
    if (!pixs || (pixGetDepth(pixs) != 32))
2140
0
        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
0
    pixGetMaxValueInRect(pixs, NULL, &maxval, NULL, NULL);
2145
0
    if (pncc) *pncc = maxval;
2146
0
    pta = ptaCreate(1);
2147
0
    ptaa = ptaaCreate(maxval + 1);
2148
0
    ptaaInitFull(ptaa, pta);
2149
0
    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
0
    pixGetDimensions(pixs, &w, &h, NULL);
2154
0
    data = pixGetData(pixs);
2155
0
    wpl = pixGetWpl(pixs);
2156
0
    for (i = 0; i < h; i++) {
2157
0
        line = data + wpl * i;
2158
0
        for (j = 0; j < w; j++) {
2159
0
            index = line[j];
2160
0
            if (index > 0)
2161
0
                ptaaAddPt(ptaa, index, j, i);
2162
0
        }
2163
0
    }
2164
2165
0
    return ptaa;
2166
0
}
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
0
{
2189
0
l_int32  w, h;
2190
0
PTA     *pta;
2191
2192
0
    if (!pixs)
2193
0
        return (PTA *)ERROR_PTR("pixs not defined", __func__, NULL);
2194
0
    pixGetDimensions(pixs, &w, &h, NULL);
2195
0
    if (x < 0 || x >= w || y < 0 || y >= h)
2196
0
        return (PTA *)ERROR_PTR("(x,y) not in pixs", __func__, NULL);
2197
0
    if (conn != 4 && conn != 8)
2198
0
        return (PTA *)ERROR_PTR("conn not 4 or 8", __func__, NULL);
2199
2200
0
    pta = ptaCreate(conn);
2201
0
    if (x > 0)
2202
0
        ptaAddPt(pta, x - 1, y);
2203
0
    if (x < w - 1)
2204
0
        ptaAddPt(pta, x + 1, y);
2205
0
    if (y > 0)
2206
0
        ptaAddPt(pta, x, y - 1);
2207
0
    if (y < h - 1)
2208
0
        ptaAddPt(pta, x, y + 1);
2209
0
    if (conn == 8) {
2210
0
        if (x > 0) {
2211
0
            if (y > 0)
2212
0
                ptaAddPt(pta, x - 1, y - 1);
2213
0
            if (y < h - 1)
2214
0
                ptaAddPt(pta, x - 1, y + 1);
2215
0
        }
2216
0
        if (x < w - 1) {
2217
0
            if (y > 0)
2218
0
                ptaAddPt(pta, x + 1, y - 1);
2219
0
            if (y < h - 1)
2220
0
                ptaAddPt(pta, x + 1, y + 1);
2221
0
        }
2222
0
    }
2223
2224
0
    return pta;
2225
0
}
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
0
{
2240
0
l_int32    i, n;
2241
0
l_float32  startx, delx, val;
2242
0
PTA       *pta;
2243
2244
0
    if (!na)
2245
0
        return (PTA *)ERROR_PTR("na not defined", __func__, NULL);
2246
2247
0
    n = numaGetCount(na);
2248
0
    pta = ptaCreate(n);
2249
0
    numaGetParameters(na, &startx, &delx);
2250
0
    for (i = 0; i < n; i++) {
2251
0
        numaGetFValue(na, i, &val);
2252
0
        ptaAddPt(pta, startx + i * delx, val);
2253
0
    }
2254
0
    return pta;
2255
0
}
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
0
{
2269
0
l_int32    i, n, nx, ny;
2270
0
l_float32  valx, valy;
2271
0
PTA       *pta;
2272
2273
0
    if (!nax || !nay)
2274
0
        return (PTA *)ERROR_PTR("nax and nay not both defined", __func__, NULL);
2275
2276
0
    nx = numaGetCount(nax);
2277
0
    ny = numaGetCount(nay);
2278
0
    n = L_MIN(nx, ny);
2279
0
    if (nx != ny)
2280
0
        L_WARNING("nx = %d does not equal ny = %d\n", __func__, nx, ny);
2281
0
    pta = ptaCreate(n);
2282
0
    for (i = 0; i < n; i++) {
2283
0
        numaGetFValue(nax, i, &valx);
2284
0
        numaGetFValue(nay, i, &valy);
2285
0
        ptaAddPt(pta, valx, valy);
2286
0
    }
2287
0
    return pta;
2288
0
}
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
0
{
2304
0
l_int32    i, n;
2305
0
l_float32  valx, valy;
2306
2307
0
    if (pnax) *pnax = NULL;
2308
0
    if (pnay) *pnay = NULL;
2309
0
    if (!pnax || !pnay)
2310
0
        return ERROR_INT("&nax and &nay not both defined", __func__, 1);
2311
0
    if (!pta)
2312
0
        return ERROR_INT("pta not defined", __func__, 1);
2313
2314
0
    n = ptaGetCount(pta);
2315
0
    *pnax = numaCreate(n);
2316
0
    *pnay = numaCreate(n);
2317
0
    for (i = 0; i < n; i++) {
2318
0
        ptaGetPt(pta, i, &valx, &valy);
2319
0
        numaAddNumber(*pnax, valx);
2320
0
        numaAddNumber(*pnay, valy);
2321
0
    }
2322
0
    return 0;
2323
0
}
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
0
{
2353
0
l_int32   i, n, w, h, x, y;
2354
0
l_uint32  rpixel, gpixel, bpixel;
2355
2356
0
    if (!pixs)
2357
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
2358
0
    if (!pta)
2359
0
        return (PIX *)ERROR_PTR("pta not defined", __func__, pixd);
2360
0
    if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
2361
0
        return (PIX *)ERROR_PTR("invalid pixd", __func__, pixd);
2362
2363
0
    if (!pixd)
2364
0
        pixd = pixConvertTo32(pixs);
2365
0
    pixGetDimensions(pixd, &w, &h, NULL);
2366
0
    composeRGBPixel(255, 0, 0, &rpixel);  /* start point */
2367
0
    composeRGBPixel(0, 255, 0, &gpixel);
2368
0
    composeRGBPixel(0, 0, 255, &bpixel);  /* end point */
2369
2370
0
    n = ptaGetCount(pta);
2371
0
    for (i = 0; i < n; i++) {
2372
0
        ptaGetIPt(pta, i, &x, &y);
2373
0
        if (x < 0 || x >= w || y < 0 || y >= h)
2374
0
            continue;
2375
0
        if (i == 0)
2376
0
            pixSetPixel(pixd, x, y, rpixel);
2377
0
        else if (i < n - 1)
2378
0
            pixSetPixel(pixd, x, y, gpixel);
2379
0
        else
2380
0
            pixSetPixel(pixd, x, y, bpixel);
2381
0
    }
2382
2383
0
    return pixd;
2384
0
}
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
0
{
2420
0
l_int32   i, n;
2421
0
l_uint32  color;
2422
0
PIXCMAP  *cmap;
2423
0
PTA      *pta;
2424
2425
0
    if (!pixs)
2426
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
2427
0
    if (!ptaa)
2428
0
        return (PIX *)ERROR_PTR("ptaa not defined", __func__, pixd);
2429
0
    if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
2430
0
        return (PIX *)ERROR_PTR("invalid pixd", __func__, pixd);
2431
0
    if (!pixp)
2432
0
        return (PIX *)ERROR_PTR("pixp not defined", __func__, pixd);
2433
2434
0
    if (!pixd)
2435
0
        pixd = pixConvertTo32(pixs);
2436
2437
        /* Use 256 random colors */
2438
0
    cmap = pixcmapCreateRandom(8, 0, 0);
2439
0
    n = ptaaGetCount(ptaa);
2440
0
    for (i = 0; i < n; i++) {
2441
0
        pixcmapGetColor32(cmap, i % 256, &color);
2442
0
        pta = ptaaGetPta(ptaa, i, L_CLONE);
2443
0
        pixDisplayPtaPattern(pixd, pixd, pta, pixp, cx, cy, color);
2444
0
        ptaDestroy(&pta);
2445
0
    }
2446
2447
0
    pixcmapDestroy(&cmap);
2448
0
    return pixd;
2449
0
}
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
0
{
2486
0
l_int32  i, n, w, h, x, y;
2487
0
PTA     *ptat;
2488
2489
0
    if (!pixs)
2490
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
2491
0
    if (!pta)
2492
0
        return (PIX *)ERROR_PTR("pta not defined", __func__, pixd);
2493
0
    if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
2494
0
        return (PIX *)ERROR_PTR("invalid pixd", __func__, pixd);
2495
0
    if (!pixp)
2496
0
        return (PIX *)ERROR_PTR("pixp not defined", __func__, pixd);
2497
2498
0
    if (!pixd)
2499
0
        pixd = pixConvertTo32(pixs);
2500
0
    pixGetDimensions(pixs, &w, &h, NULL);
2501
0
    ptat = ptaReplicatePattern(pta, pixp, NULL, cx, cy, w, h);
2502
2503
0
    n = ptaGetCount(ptat);
2504
0
    for (i = 0; i < n; i++) {
2505
0
        ptaGetIPt(ptat, i, &x, &y);
2506
0
        if (x < 0 || x >= w || y < 0 || y >= h)
2507
0
            continue;
2508
0
        pixSetPixel(pixd, x, y, color);
2509
0
    }
2510
2511
0
    ptaDestroy(&ptat);
2512
0
    return pixd;
2513
0
}
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
0
{
2545
0
l_int32  i, j, n, np, x, y, xp, yp, xf, yf;
2546
0
PTA     *ptat, *ptad;
2547
2548
0
    if (!ptas)
2549
0
        return (PTA *)ERROR_PTR("ptas not defined", __func__, NULL);
2550
0
    if (!pixp && !ptap)
2551
0
        return (PTA *)ERROR_PTR("no pattern is defined", __func__, NULL);
2552
0
    if (pixp && ptap)
2553
0
        L_WARNING("pixp and ptap defined; using ptap\n", __func__);
2554
2555
0
    n = ptaGetCount(ptas);
2556
0
    ptad = ptaCreate(n);
2557
0
    if (ptap)
2558
0
        ptat = ptaClone(ptap);
2559
0
    else
2560
0
        ptat = ptaGetPixelsFromPix(pixp, NULL);
2561
0
    np = ptaGetCount(ptat);
2562
0
    for (i = 0; i < n; i++) {
2563
0
        ptaGetIPt(ptas, i, &x, &y);
2564
0
        for (j = 0; j < np; j++) {
2565
0
            ptaGetIPt(ptat, j, &xp, &yp);
2566
0
            xf = x - cx + xp;
2567
0
            yf = y - cy + yp;
2568
0
            if (xf >= 0 && xf < w && yf >= 0 && yf < h)
2569
0
                ptaAddPt(ptad, xf, yf);
2570
0
        }
2571
0
    }
2572
2573
0
    ptaDestroy(&ptat);
2574
0
    return ptad;
2575
0
}
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
0
{
2590
0
l_int32    i, j, w, h, npta, npt, x, y, rv, gv, bv;
2591
0
l_uint32  *pixela;
2592
0
NUMA      *na1, *na2, *na3;
2593
0
PIX       *pixd;
2594
0
PTA       *pta;
2595
2596
0
    if (!pixs)
2597
0
        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
2598
0
    if (!ptaa)
2599
0
        return (PIX *)ERROR_PTR("ptaa not defined", __func__, NULL);
2600
0
    npta = ptaaGetCount(ptaa);
2601
0
    if (npta == 0)
2602
0
        return (PIX *)ERROR_PTR("no pta", __func__, NULL);
2603
2604
0
    if ((pixd = pixConvertTo32(pixs)) == NULL)
2605
0
        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
2606
0
    pixGetDimensions(pixd, &w, &h, NULL);
2607
2608
        /* Make a colormap for the paths */
2609
0
    if ((pixela = (l_uint32 *)LEPT_CALLOC(npta, sizeof(l_uint32))) == NULL) {
2610
0
        pixDestroy(&pixd);
2611
0
        return (PIX *)ERROR_PTR("calloc fail for pixela", __func__, NULL);
2612
0
    }
2613
0
    na1 = numaPseudorandomSequence(256, 14657);
2614
0
    na2 = numaPseudorandomSequence(256, 34631);
2615
0
    na3 = numaPseudorandomSequence(256, 54617);
2616
0
    for (i = 0; i < npta; i++) {
2617
0
        numaGetIValue(na1, i % 256, &rv);
2618
0
        numaGetIValue(na2, i % 256, &gv);
2619
0
        numaGetIValue(na3, i % 256, &bv);
2620
0
        composeRGBPixel(rv, gv, bv, &pixela[i]);
2621
0
    }
2622
0
    numaDestroy(&na1);
2623
0
    numaDestroy(&na2);
2624
0
    numaDestroy(&na3);
2625
2626
0
    for (i = 0; i < npta; i++) {
2627
0
        pta = ptaaGetPta(ptaa, i, L_CLONE);
2628
0
        npt = ptaGetCount(pta);
2629
0
        for (j = 0; j < npt; j++) {
2630
0
            ptaGetIPt(pta, j, &x, &y);
2631
0
            if (x < 0 || x >= w || y < 0 || y >= h)
2632
0
                continue;
2633
0
            pixSetPixel(pixd, x, y, pixela[i]);
2634
0
        }
2635
0
        ptaDestroy(&pta);
2636
0
    }
2637
2638
0
    LEPT_FREE(pixela);
2639
0
    return pixd;
2640
0
}