/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 | } |