Mercurial > hgrepos > Python2 > PyMuPDF

comparison mupdf-source/thirdparty/leptonica/src/affine.c @ 2:b50eed0cc0ef upstream

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ADD: MuPDF v1.26.7: the MuPDF source as downloaded by a default build of PyMuPDF 1.26.4. The directory name has changed: no version number in the expanded directory now.
author Franz Glasner <fzglas.hg@dom66.de>
date Mon, 15 Sep 2025 11:43:07 +0200
parents
children
comparison
equal deleted inserted replaced
1:1d09e1dec1d9 2:b50eed0cc0ef
1 /*====================================================================*
2 - Copyright (C) 2001 Leptonica. All rights reserved.
3 -
4 - Redistribution and use in source and binary forms, with or without
5 - modification, are permitted provided that the following conditions
6 - are met:
7 - 1. Redistributions of source code must retain the above copyright
8 - notice, this list of conditions and the following disclaimer.
9 - 2. Redistributions in binary form must reproduce the above
10 - copyright notice, this list of conditions and the following
11 - disclaimer in the documentation and/or other materials
12 - provided with the distribution.
13 -
14 - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
15 - ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
16 - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
17 - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ANY
18 - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
19 - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
20 - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
21 - PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
22 - OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
23 - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 *====================================================================*/
26
27
28 /*!
29 * \file affine.c
30 * <pre>
31 *
32 * Affine (3 pt) image transformation using a sampled
33 * (to nearest integer) transform on each dest point
34 * PIX *pixAffineSampledPta()
35 * PIX *pixAffineSampled()
36 *
37 * Affine (3 pt) image transformation using interpolation
38 * (or area mapping) for anti-aliasing images that are
39 * 2, 4, or 8 bpp gray, or colormapped, or 32 bpp RGB
40 * PIX *pixAffinePta()
41 * PIX *pixAffine()
42 * PIX *pixAffinePtaColor()
43 * PIX *pixAffineColor()
44 * PIX *pixAffinePtaGray()
45 * PIX *pixAffineGray()
46 *
47 * Affine transform including alpha (blend) component
48 * PIX *pixAffinePtaWithAlpha()
49 *
50 * Affine coordinate transformation
51 * l_int32 getAffineXformCoeffs()
52 * l_int32 affineInvertXform()
53 * l_int32 affineXformSampledPt()
54 * l_int32 affineXformPt()
55 *
56 * Interpolation helper functions
57 * l_int32 linearInterpolatePixelGray()
58 * l_int32 linearInterpolatePixelColor()
59 *
60 * Gauss-jordan linear equation solver
61 * l_int32 gaussjordan()
62 *
63 * Affine image transformation using a sequence of
64 * shear/scale/translation operations
65 * PIX *pixAffineSequential()
66 *
67 * One can define a coordinate space by the location of the origin,
68 * the orientation of x and y axes, and the unit scaling along
69 * each axis. An affine transform is a general linear
70 * transformation from one coordinate space to another.
71 *
72 * For the general case, we can define the affine transform using
73 * two sets of three (noncollinear) points in a plane. One set
74 * corresponds to the input (src) coordinate space; the other to the
75 * transformed (dest) coordinate space. Each point in the
76 * src corresponds to one of the points in the dest. With two
77 * sets of three points, we get a set of 6 equations in 6 unknowns
78 * that specifies the mapping between the coordinate spaces.
79 * The interface here allows you to specify either the corresponding
80 * sets of 3 points, or the transform itself (as a vector of 6
81 * coefficients).
82 *
83 * Given the transform as a vector of 6 coefficients, we can compute
84 * both a a pointwise affine coordinate transformation and an
85 * affine image transformation.
86 *
87 * To compute the coordinate transform, we need the coordinate
88 * value (x',y') in the transformed space for any point (x,y)
89 * in the original space. To derive this transform from the
90 * three corresponding points, it is convenient to express the affine
91 * coordinate transformation using an LU decomposition of
92 * a set of six linear equations that express the six coordinates
93 * of the three points in the transformed space as a function of
94 * the six coordinates in the original space. Once we have
95 * this transform matrix , we can transform an image by
96 * finding, for each destination pixel, the pixel (or pixels)
97 * in the source that give rise to it.
98 *
99 * This 'pointwise' transformation can be done either by sampling
100 * and picking a single pixel in the src to replicate into the dest,
101 * or by interpolating (or averaging) over four src pixels to
102 * determine the value of the dest pixel. The first method is
103 * implemented by pixAffineSampled() and the second method by
104 * pixAffine(). The interpolated method can only be used for
105 * images with more than 1 bpp, but for these, the image quality
106 * is significantly better than the sampled method, due to
107 * the 'antialiasing' effect of weighting the src pixels.
108 *
109 * Interpolation works well when there is relatively little scaling,
110 * or if there is image expansion in general. However, if there
111 * is significant image reduction, one should apply a low-pass
112 * filter before subsampling to avoid aliasing the high frequencies.
113 *
114 * A typical application might be to align two images, which
115 * may be scaled, rotated and translated versions of each other.
116 * Through some pre-processing, three corresponding points are
117 * located in each of the two images. One of the images is
118 * then to be (affine) transformed to align with the other.
119 * As mentioned, the standard way to do this is to use three
120 * sets of points, compute the 6 transformation coefficients
121 * from these points that describe the linear transformation,
122 *
123 * x' = ax + by + c
124 * y' = dx + ey + f
125 *
126 * and use this in a pointwise manner to transform the image.
127 *
128 * N.B. Be sure to see the comment in getAffineXformCoeffs(),
129 * regarding using the inverse of the affine transform for points
130 * to transform images.
131 *
132 * There is another way to do this transformation; namely,
133 * by doing a sequence of simple affine transforms, without
134 * computing directly the affine coordinate transformation.
135 * We have at our disposal (1) translations (using rasterop),
136 * (2) horizontal and vertical shear about any horizontal and vertical
137 * line, respectively, and (3) non-isotropic scaling by two
138 * arbitrary x and y scaling factors. We also have rotation
139 * about an arbitrary point, but this is equivalent to a set
140 * of three shears so we do not need to use it.
141 *
142 * Why might we do this? For binary images, it is usually
143 * more efficient to do such transformations by a sequence
144 * of word parallel operations. Shear and translation can be
145 * done in-place and word parallel; arbitrary scaling is
146 * mostly pixel-wise.
147 *
148 * Suppose that we are transforming image 1 to correspond to image 2.
149 * We have a set of three points, describing the coordinate space
150 * embedded in image 1, and we need to transform image 1 until
151 * those three points exactly correspond to the new coordinate space
152 * defined by the second set of three points. In our image
153 * matching application, the latter set of three points was
154 * found to be the corresponding points in image 2.
155 *
156 * The most elegant way I can think of to do such a sequential
157 * implementation is to imagine that we're going to transform
158 * BOTH images until they're aligned. (We don't really want
159 * to transform both, because in fact we may only have one image
160 * that is undergoing a general affine transformation.)
161 *
162 * Choose the 3 corresponding points as follows:
163 * ~ The 1st point is an origin
164 * ~ The 2nd point gives the orientation and scaling of the
165 * "x" axis with respect to the origin
166 * ~ The 3rd point does likewise for the "y" axis.
167 * These "axes" must not be collinear; otherwise they are
168 * arbitrary (although some strange things will happen if
169 * the handedness sweeping through the minimum angle between
170 * the axes is opposite).
171 *
172 * An important constraint is that we have shear operations
173 * about an arbitrary horizontal or vertical line, but always
174 * parallel to the x or y axis. If we continue to pretend that
175 * we have an unprimed coordinate space embedded in image 1 and
176 * a primed coordinate space embedded in image 2, we imagine
177 * (a) transforming image 1 by horizontal and vertical shears about
178 * point 1 to align points 3 and 2 along the y and x axes,
179 * respectively, and (b) transforming image 2 by horizontal and
180 * vertical shears about point 1' to align points 3' and 2' along
181 * the y and x axes. Then we scale image 1 so that the distances
182 * from 1 to 2 and from 1 to 3 are equal to the distances in
183 * image 2 from 1' to 2' and from 1' to 3'. This scaling operation
184 * leaves the true image origin, at (0,0) invariant, and will in
185 * general translate point 1. The original points 1 and 1' will
186 * typically not coincide in any event, so we must translate
187 * the origin of image 1, at its current point 1, to the origin
188 * of image 2 at 1'. The images should now be aligned. But
189 * because we never really transformed image 2 (and image 2 may
190 * not even exist), we now perform on image 1 the reverse of
191 * the shear transforms that we imagined doing on image 2;
192 * namely, the negative vertical shear followed by the negative
193 * horizontal shear. Image 1 should now have its transformed
194 * unprimed coordinates aligned with the original primed
195 * coordinates. In all this, it is only necessary to keep track
196 * of the shear angles and translations of points during the shears.
197 * What has been accomplished is a general affine transformation
198 * on image 1.
199 *
200 * Having described all this, if you are going to use an
201 * affine transformation in an application, this is what you
202 * need to know:
203 *
204 * (1) You should NEVER use the sequential method, because
205 * the image quality for 1 bpp text is much poorer
206 * (even though it is about 2x faster than the pointwise sampled
207 * method), and for images with depth greater than 1, it is
208 * nearly 20x slower than the pointwise sampled method
209 * and over 10x slower than the pointwise interpolated method!
210 * The sequential method is given here for purely
211 * pedagogical reasons.
212 *
213 * (2) For 1 bpp images, use the pointwise sampled function
214 * pixAffineSampled(). For all other images, the best
215 * quality results result from using the pointwise
216 * interpolated function pixAffinePta() or pixAffine();
217 * the cost is less than a doubling of the computation time
218 * with respect to the sampled function. If you use
219 * interpolation on colormapped images, the colormap will
220 * be removed, resulting in either a grayscale or color
221 * image, depending on the values in the colormap.
222 * If you want to retain the colormap, use pixAffineSampled().
223 *
224 * Typical relative timing of pointwise transforms (sampled = 1.0):
225 * 8 bpp: sampled 1.0
226 * interpolated 1.6
227 * 32 bpp: sampled 1.0
228 * interpolated 1.8
229 * Additionally, the computation time/pixel is nearly the same
230 * for 8 bpp and 32 bpp, for both sampled and interpolated.
231 * </pre>
232 */
233
234 #ifdef HAVE_CONFIG_H
235 #include <config_auto.h>
236 #endif /* HAVE_CONFIG_H */
237
238 #include <string.h>
239 #include <math.h>
240 #include "allheaders.h"
241
242 extern l_float32 AlphaMaskBorderVals[2];
243
244 #ifndef NO_CONSOLE_IO
245 #define DEBUG 0
246 #endif /* ~NO_CONSOLE_IO */
247
248 /*-------------------------------------------------------------*
249 * Sampled affine image transformation *
250 *-------------------------------------------------------------*/
251 /*!
252 * \brief pixAffineSampledPta()
253 *
254 * \param[in] pixs all depths
255 * \param[in] ptad 3 pts of final coordinate space
256 * \param[in] ptas 3 pts of initial coordinate space
257 * \param[in] incolor L_BRING_IN_WHITE, L_BRING_IN_BLACK
258 * \return pixd, or NULL on error
259 *
260 * <pre>
261 * Notes:
262 * (1) Brings in either black or white pixels from the boundary.
263 * (2) Retains colormap, which you can do for a sampled transform..
264 * (3) The 3 points must not be collinear.
265 * (4) The order of the 3 points is arbitrary; however, to compare
266 * with the sequential transform they must be in these locations
267 * and in this order: origin, x-axis, y-axis.
268 * (5) For 1 bpp images, this has much better quality results
269 * than pixAffineSequential(), particularly for text.
270 * It is about 3x slower, but does not require additional
271 * border pixels. The poor quality of pixAffineSequential()
272 * is due to repeated quantized transforms. It is strongly
273 * recommended that pixAffineSampled() be used for 1 bpp images.
274 * (6) For 8 or 32 bpp, much better quality is obtained by the
275 * somewhat slower pixAffinePta(). See that function
276 * for relative timings between sampled and interpolated.
277 * (7) To repeat, use of the sequential transform,
278 * pixAffineSequential(), for any images, is discouraged.
279 * </pre>
280 */
281 PIX *
282 pixAffineSampledPta(PIX *pixs,
283 PTA *ptad,
284 PTA *ptas,
285 l_int32 incolor)
286 {
287 l_float32 *vc;
288 PIX *pixd;
289
290 if (!pixs)
291 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
292 if (!ptas)
293 return (PIX *)ERROR_PTR("ptas not defined", __func__, NULL);
294 if (!ptad)
295 return (PIX *)ERROR_PTR("ptad not defined", __func__, NULL);
296 if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
297 return (PIX *)ERROR_PTR("invalid incolor", __func__, NULL);
298 if (ptaGetCount(ptas) != 3)
299 return (PIX *)ERROR_PTR("ptas count not 3", __func__, NULL);
300 if (ptaGetCount(ptad) != 3)
301 return (PIX *)ERROR_PTR("ptad count not 3", __func__, NULL);
302
303 /* Get backwards transform from dest to src, and apply it */
304 getAffineXformCoeffs(ptad, ptas, &vc);
305 pixd = pixAffineSampled(pixs, vc, incolor);
306 LEPT_FREE(vc);
307
308 return pixd;
309 }
310
311
312 /*!
313 * \brief pixAffineSampled()
314 *
315 * \param[in] pixs all depths
316 * \param[in] vc vector of 6 coefficients for affine transformation
317 * \param[in] incolor L_BRING_IN_WHITE, L_BRING_IN_BLACK
318 * \return pixd, or NULL on error
319 *
320 * <pre>
321 * Notes:
322 * (1) Brings in either black or white pixels from the boundary.
323 * (2) Retains colormap, which you can do for a sampled transform..
324 * (3) For 8 or 32 bpp, much better quality is obtained by the
325 * somewhat slower pixAffine(). See that function
326 * for relative timings between sampled and interpolated.
327 * </pre>
328 */
329 PIX *
330 pixAffineSampled(PIX *pixs,
331 l_float32 *vc,
332 l_int32 incolor)
333 {
334 l_int32 i, j, w, h, d, x, y, wpls, wpld, color, cmapindex;
335 l_uint32 val;
336 l_uint32 *datas, *datad, *lines, *lined;
337 PIX *pixd;
338 PIXCMAP *cmap;
339
340 if (!pixs)
341 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
342 if (!vc)
343 return (PIX *)ERROR_PTR("vc not defined", __func__, NULL);
344 if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
345 return (PIX *)ERROR_PTR("invalid incolor", __func__, NULL);
346 pixGetDimensions(pixs, &w, &h, &d);
347 if (d != 1 && d != 2 && d != 4 && d != 8 && d != 32)
348 return (PIX *)ERROR_PTR("depth not 1, 2, 4, 8 or 16", __func__, NULL);
349
350 /* Init all dest pixels to color to be brought in from outside */
351 pixd = pixCreateTemplate(pixs);
352 if ((cmap = pixGetColormap(pixs)) != NULL) {
353 if (incolor == L_BRING_IN_WHITE)
354 color = 1;
355 else
356 color = 0;
357 pixcmapAddBlackOrWhite(cmap, color, &cmapindex);
358 pixSetAllArbitrary(pixd, cmapindex);
359 } else {
360 if ((d == 1 && incolor == L_BRING_IN_WHITE) ||
361 (d > 1 && incolor == L_BRING_IN_BLACK)) {
362 pixClearAll(pixd);
363 } else {
364 pixSetAll(pixd);
365 }
366 }
367
368 /* Scan over the dest pixels */
369 datas = pixGetData(pixs);
370 wpls = pixGetWpl(pixs);
371 datad = pixGetData(pixd);
372 wpld = pixGetWpl(pixd);
373 for (i = 0; i < h; i++) {
374 lined = datad + i * wpld;
375 for (j = 0; j < w; j++) {
376 affineXformSampledPt(vc, j, i, &x, &y);
377 if (x < 0 || y < 0 || x >=w || y >= h)
378 continue;
379 lines = datas + y * wpls;
380 if (d == 1) {
381 val = GET_DATA_BIT(lines, x);
382 SET_DATA_BIT_VAL(lined, j, val);
383 } else if (d == 8) {
384 val = GET_DATA_BYTE(lines, x);
385 SET_DATA_BYTE(lined, j, val);
386 } else if (d == 32) {
387 lined[j] = lines[x];
388 } else if (d == 2) {
389 val = GET_DATA_DIBIT(lines, x);
390 SET_DATA_DIBIT(lined, j, val);
391 } else if (d == 4) {
392 val = GET_DATA_QBIT(lines, x);
393 SET_DATA_QBIT(lined, j, val);
394 }
395 }
396 }
397
398 return pixd;
399 }
400
401
402 /*---------------------------------------------------------------------*
403 * Interpolated affine image transformation *
404 *---------------------------------------------------------------------*/
405 /*!
406 * \brief pixAffinePta()
407 *
408 * \param[in] pixs all depths; colormap ok
409 * \param[in] ptad 3 pts of final coordinate space
410 * \param[in] ptas 3 pts of initial coordinate space
411 * \param[in] incolor L_BRING_IN_WHITE, L_BRING_IN_BLACK
412 * \return pixd, or NULL on error
413 *
414 * <pre>
415 * Notes:
416 * (1) Brings in either black or white pixels from the boundary
417 * (2) Removes any existing colormap, if necessary, before transforming
418 * </pre>
419 */
420 PIX *
421 pixAffinePta(PIX *pixs,
422 PTA *ptad,
423 PTA *ptas,
424 l_int32 incolor)
425 {
426 l_int32 d;
427 l_uint32 colorval;
428 PIX *pixt1, *pixt2, *pixd;
429
430 if (!pixs)
431 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
432 if (!ptas)
433 return (PIX *)ERROR_PTR("ptas not defined", __func__, NULL);
434 if (!ptad)
435 return (PIX *)ERROR_PTR("ptad not defined", __func__, NULL);
436 if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
437 return (PIX *)ERROR_PTR("invalid incolor", __func__, NULL);
438 if (ptaGetCount(ptas) != 3)
439 return (PIX *)ERROR_PTR("ptas count not 3", __func__, NULL);
440 if (ptaGetCount(ptad) != 3)
441 return (PIX *)ERROR_PTR("ptad count not 3", __func__, NULL);
442
443 if (pixGetDepth(pixs) == 1)
444 return pixAffineSampledPta(pixs, ptad, ptas, incolor);
445
446 /* Remove cmap if it exists, and unpack to 8 bpp if necessary */
447 pixt1 = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
448 d = pixGetDepth(pixt1);
449 if (d < 8)
450 pixt2 = pixConvertTo8(pixt1, FALSE);
451 else
452 pixt2 = pixClone(pixt1);
453 d = pixGetDepth(pixt2);
454
455 /* Compute actual color to bring in from edges */
456 colorval = 0;
457 if (incolor == L_BRING_IN_WHITE) {
458 if (d == 8)
459 colorval = 255;
460 else /* d == 32 */
461 colorval = 0xffffff00;
462 }
463
464 if (d == 8)
465 pixd = pixAffinePtaGray(pixt2, ptad, ptas, colorval);
466 else /* d == 32 */
467 pixd = pixAffinePtaColor(pixt2, ptad, ptas, colorval);
468 pixDestroy(&pixt1);
469 pixDestroy(&pixt2);
470 return pixd;
471 }
472
473
474 /*!
475 * \brief pixAffine()
476 *
477 * \param[in] pixs all depths; colormap ok
478 * \param[in] vc vector of 6 coefficients for affine transformation
479 * \param[in] incolor L_BRING_IN_WHITE, L_BRING_IN_BLACK
480 * \return pixd, or NULL on error
481 *
482 * <pre>
483 * Notes:
484 * (1) Brings in either black or white pixels from the boundary
485 * (2) Removes any existing colormap, if necessary, before transforming
486 * </pre>
487 */
488 PIX *
489 pixAffine(PIX *pixs,
490 l_float32 *vc,
491 l_int32 incolor)
492 {
493 l_int32 d;
494 l_uint32 colorval;
495 PIX *pixt1, *pixt2, *pixd;
496
497 if (!pixs)
498 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
499 if (!vc)
500 return (PIX *)ERROR_PTR("vc not defined", __func__, NULL);
501
502 if (pixGetDepth(pixs) == 1)
503 return pixAffineSampled(pixs, vc, incolor);
504
505 /* Remove cmap if it exists, and unpack to 8 bpp if necessary */
506 pixt1 = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
507 d = pixGetDepth(pixt1);
508 if (d < 8)
509 pixt2 = pixConvertTo8(pixt1, FALSE);
510 else
511 pixt2 = pixClone(pixt1);
512 d = pixGetDepth(pixt2);
513
514 /* Compute actual color to bring in from edges */
515 colorval = 0;
516 if (incolor == L_BRING_IN_WHITE) {
517 if (d == 8)
518 colorval = 255;
519 else /* d == 32 */
520 colorval = 0xffffff00;
521 }
522
523 if (d == 8)
524 pixd = pixAffineGray(pixt2, vc, colorval);
525 else /* d == 32 */
526 pixd = pixAffineColor(pixt2, vc, colorval);
527 pixDestroy(&pixt1);
528 pixDestroy(&pixt2);
529 return pixd;
530 }
531
532
533 /*!
534 * \brief pixAffinePtaColor()
535 *
536 * \param[in] pixs 32 bpp
537 * \param[in] ptad 3 pts of final coordinate space
538 * \param[in] ptas 3 pts of initial coordinate space
539 * \param[in] colorval e.g.: 0 to bring in BLACK, 0xffffff00 for WHITE
540 * \return pixd, or NULL on error
541 */
542 PIX *
543 pixAffinePtaColor(PIX *pixs,
544 PTA *ptad,
545 PTA *ptas,
546 l_uint32 colorval)
547 {
548 l_float32 *vc;
549 PIX *pixd;
550
551 if (!pixs)
552 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
553 if (!ptas)
554 return (PIX *)ERROR_PTR("ptas not defined", __func__, NULL);
555 if (!ptad)
556 return (PIX *)ERROR_PTR("ptad not defined", __func__, NULL);
557 if (pixGetDepth(pixs) != 32)
558 return (PIX *)ERROR_PTR("pixs must be 32 bpp", __func__, NULL);
559 if (ptaGetCount(ptas) != 3)
560 return (PIX *)ERROR_PTR("ptas count not 3", __func__, NULL);
561 if (ptaGetCount(ptad) != 3)
562 return (PIX *)ERROR_PTR("ptad count not 3", __func__, NULL);
563
564 /* Get backwards transform from dest to src, and apply it */
565 getAffineXformCoeffs(ptad, ptas, &vc);
566 pixd = pixAffineColor(pixs, vc, colorval);
567 LEPT_FREE(vc);
568
569 return pixd;
570 }
571
572
573 /*!
574 * \brief pixAffineColor()
575 *
576 * \param[in] pixs 32 bpp
577 * \param[in] vc vector of 6 coefficients for affine transformation
578 * \param[in] colorval e.g.: 0 to bring in BLACK, 0xffffff00 for WHITE
579 * \return pixd, or NULL on error
580 */
581 PIX *
582 pixAffineColor(PIX *pixs,
583 l_float32 *vc,
584 l_uint32 colorval)
585 {
586 l_int32 i, j, w, h, d, wpls, wpld;
587 l_uint32 val;
588 l_uint32 *datas, *datad, *lined;
589 l_float32 x, y;
590 PIX *pix1, *pix2, *pixd;
591
592 if (!pixs)
593 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
594 pixGetDimensions(pixs, &w, &h, &d);
595 if (d != 32)
596 return (PIX *)ERROR_PTR("pixs must be 32 bpp", __func__, NULL);
597 if (!vc)
598 return (PIX *)ERROR_PTR("vc not defined", __func__, NULL);
599
600 datas = pixGetData(pixs);
601 wpls = pixGetWpl(pixs);
602 pixd = pixCreateTemplate(pixs);
603 pixSetAllArbitrary(pixd, colorval);
604 datad = pixGetData(pixd);
605 wpld = pixGetWpl(pixd);
606
607 /* Iterate over destination pixels */
608 for (i = 0; i < h; i++) {
609 lined = datad + i * wpld;
610 for (j = 0; j < w; j++) {
611 /* Compute float src pixel location corresponding to (i,j) */
612 affineXformPt(vc, j, i, &x, &y);
613 linearInterpolatePixelColor(datas, wpls, w, h, x, y, colorval,
614 &val);
615 *(lined + j) = val;
616 }
617 }
618
619 /* If rgba, transform the pixs alpha channel and insert in pixd */
620 if (pixGetSpp(pixs) == 4) {
621 pix1 = pixGetRGBComponent(pixs, L_ALPHA_CHANNEL);
622 pix2 = pixAffineGray(pix1, vc, 255); /* bring in opaque */
623 pixSetRGBComponent(pixd, pix2, L_ALPHA_CHANNEL);
624 pixDestroy(&pix1);
625 pixDestroy(&pix2);
626 }
627
628 return pixd;
629 }
630
631
632 /*!
633 * \brief pixAffinePtaGray()
634 *
635 * \param[in] pixs 8 bpp
636 * \param[in] ptad 3 pts of final coordinate space
637 * \param[in] ptas 3 pts of initial coordinate space
638 * \param[in] grayval e.g.: 0 to bring in BLACK, 255 for WHITE
639 * \return pixd, or NULL on error
640 */
641 PIX *
642 pixAffinePtaGray(PIX *pixs,
643 PTA *ptad,
644 PTA *ptas,
645 l_uint8 grayval)
646 {
647 l_float32 *vc;
648 PIX *pixd;
649
650 if (!pixs)
651 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
652 if (!ptas)
653 return (PIX *)ERROR_PTR("ptas not defined", __func__, NULL);
654 if (!ptad)
655 return (PIX *)ERROR_PTR("ptad not defined", __func__, NULL);
656 if (pixGetDepth(pixs) != 8)
657 return (PIX *)ERROR_PTR("pixs must be 8 bpp", __func__, NULL);
658 if (ptaGetCount(ptas) != 3)
659 return (PIX *)ERROR_PTR("ptas count not 3", __func__, NULL);
660 if (ptaGetCount(ptad) != 3)
661 return (PIX *)ERROR_PTR("ptad count not 3", __func__, NULL);
662
663 /* Get backwards transform from dest to src, and apply it */
664 getAffineXformCoeffs(ptad, ptas, &vc);
665 pixd = pixAffineGray(pixs, vc, grayval);
666 LEPT_FREE(vc);
667
668 return pixd;
669 }
670
671
672
673 /*!
674 * \brief pixAffineGray()
675 *
676 * \param[in] pixs 8 bpp
677 * \param[in] vc vector of 6 coefficients for affine transformation
678 * \param[in] grayval e.g.: 0 to bring in BLACK, 255 for WHITE
679 * \return pixd, or NULL on error
680 */
681 PIX *
682 pixAffineGray(PIX *pixs,
683 l_float32 *vc,
684 l_uint8 grayval)
685 {
686 l_int32 i, j, w, h, wpls, wpld, val;
687 l_uint32 *datas, *datad, *lined;
688 l_float32 x, y;
689 PIX *pixd;
690
691 if (!pixs)
692 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
693 pixGetDimensions(pixs, &w, &h, NULL);
694 if (pixGetDepth(pixs) != 8)
695 return (PIX *)ERROR_PTR("pixs must be 8 bpp", __func__, NULL);
696 if (!vc)
697 return (PIX *)ERROR_PTR("vc not defined", __func__, NULL);
698
699 datas = pixGetData(pixs);
700 wpls = pixGetWpl(pixs);
701 pixd = pixCreateTemplate(pixs);
702 pixSetAllArbitrary(pixd, grayval);
703 datad = pixGetData(pixd);
704 wpld = pixGetWpl(pixd);
705
706 /* Iterate over destination pixels */
707 for (i = 0; i < h; i++) {
708 lined = datad + i * wpld;
709 for (j = 0; j < w; j++) {
710 /* Compute float src pixel location corresponding to (i,j) */
711 affineXformPt(vc, j, i, &x, &y);
712 linearInterpolatePixelGray(datas, wpls, w, h, x, y, grayval, &val);
713 SET_DATA_BYTE(lined, j, val);
714 }
715 }
716
717 return pixd;
718 }
719
720
721 /*---------------------------------------------------------------------------*
722 * Affine transform including alpha (blend) component *
723 *---------------------------------------------------------------------------*/
724 /*!
725 * \brief pixAffinePtaWithAlpha()
726 *
727 * \param[in] pixs 32 bpp rgb
728 * \param[in] ptad 3 pts of final coordinate space
729 * \param[in] ptas 3 pts of initial coordinate space
730 * \param[in] pixg [optional] 8 bpp, can be null
731 * \param[in] fract between 0.0 and 1.0, with 0.0 fully transparent
732 * and 1.0 fully opaque
733 * \param[in] border of pixels added to capture transformed source pixels
734 * \return pixd, or NULL on error
735 *
736 * <pre>
737 * Notes:
738 * (1) The alpha channel is transformed separately from pixs,
739 * and aligns with it, being fully transparent outside the
740 * boundary of the transformed pixs. For pixels that are fully
741 * transparent, a blending function like pixBlendWithGrayMask()
742 * will give zero weight to corresponding pixels in pixs.
743 * (2) If pixg is NULL, it is generated as an alpha layer that is
744 * partially opaque, using %fract. Otherwise, it is cropped
745 * to pixs if required and %fract is ignored. The alpha channel
746 * in pixs is never used.
747 * (3) Colormaps are removed.
748 * (4) When pixs is transformed, it doesn't matter what color is brought
749 * in because the alpha channel will be transparent (0) there.
750 * (5) To avoid losing source pixels in the destination, it may be
751 * necessary to add a border to the source pix before doing
752 * the affine transformation. This can be any non-negative number.
753 * (6) The input %ptad and %ptas are in a coordinate space before
754 * the border is added. Internally, we compensate for this
755 * before doing the affine transform on the image after the border
756 * is added.
757 * (7) The default setting for the border values in the alpha channel
758 * is 0 (transparent) for the outermost ring of pixels and
759 * (0.5 * fract * 255) for the second ring. When blended over
760 * a second image, this
761 * (a) shrinks the visible image to make a clean overlap edge
762 * with an image below, and
763 * (b) softens the edges by weakening the aliasing there.
764 * Use l_setAlphaMaskBorder() to change these values.
765 * </pre>
766 */
767 PIX *
768 pixAffinePtaWithAlpha(PIX *pixs,
769 PTA *ptad,
770 PTA *ptas,
771 PIX *pixg,
772 l_float32 fract,
773 l_int32 border)
774 {
775 l_int32 ws, hs, d;
776 PIX *pixd, *pixb1, *pixb2, *pixg2, *pixga;
777 PTA *ptad2, *ptas2;
778
779 if (!pixs)
780 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
781 pixGetDimensions(pixs, &ws, &hs, &d);
782 if (d != 32 && pixGetColormap(pixs) == NULL)
783 return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", __func__, NULL);
784 if (pixg && pixGetDepth(pixg) != 8) {
785 L_WARNING("pixg not 8 bpp; using 'fract' transparent alpha\n",
786 __func__);
787 pixg = NULL;
788 }
789 if (!pixg && (fract < 0.0 || fract > 1.0)) {
790 L_WARNING("invalid fract; using 1.0 (fully transparent)\n", __func__);
791 fract = 1.0;
792 }
793 if (!pixg && fract == 0.0)
794 L_WARNING("fully opaque alpha; image will not be blended\n", __func__);
795 if (!ptad)
796 return (PIX *)ERROR_PTR("ptad not defined", __func__, NULL);
797 if (!ptas)
798 return (PIX *)ERROR_PTR("ptas not defined", __func__, NULL);
799
800 /* Add border; the color doesn't matter */
801 pixb1 = pixAddBorder(pixs, border, 0);
802
803 /* Transform the ptr arrays to work on the bordered image */
804 ptad2 = ptaTransform(ptad, border, border, 1.0, 1.0);
805 ptas2 = ptaTransform(ptas, border, border, 1.0, 1.0);
806
807 /* Do separate affine transform of rgb channels of pixs and of pixg */
808 pixd = pixAffinePtaColor(pixb1, ptad2, ptas2, 0);
809 if (!pixg) {
810 pixg2 = pixCreate(ws, hs, 8);
811 if (fract == 1.0)
812 pixSetAll(pixg2);
813 else
814 pixSetAllArbitrary(pixg2, (l_int32)(255.0 * fract));
815 } else {
816 pixg2 = pixResizeToMatch(pixg, NULL, ws, hs);
817 }
818 if (ws > 10 && hs > 10) { /* see note 7 */
819 pixSetBorderRingVal(pixg2, 1,
820 (l_int32)(255.0 * fract * AlphaMaskBorderVals[0]));
821 pixSetBorderRingVal(pixg2, 2,
822 (l_int32)(255.0 * fract * AlphaMaskBorderVals[1]));
823
824 }
825 pixb2 = pixAddBorder(pixg2, border, 0); /* must be black border */
826 pixga = pixAffinePtaGray(pixb2, ptad2, ptas2, 0);
827 pixSetRGBComponent(pixd, pixga, L_ALPHA_CHANNEL);
828 pixSetSpp(pixd, 4);
829
830 pixDestroy(&pixg2);
831 pixDestroy(&pixb1);
832 pixDestroy(&pixb2);
833 pixDestroy(&pixga);
834 ptaDestroy(&ptad2);
835 ptaDestroy(&ptas2);
836 return pixd;
837 }
838
839
840 /*-------------------------------------------------------------*
841 * Affine coordinate transformation *
842 *-------------------------------------------------------------*/
843 /*!
844 * \brief getAffineXformCoeffs()
845 *
846 * \param[in] ptas source 3 points; unprimed
847 * \param[in] ptad transformed 3 points; primed
848 * \param[out] pvc vector of coefficients of transform
849 * \return 0 if OK; 1 on error
850 *
851 * <pre>
852 * We have a set of six equations, describing the affine
853 * transformation that takes 3 points ptas into 3 other
854 * points ptad. These equations are:
855 *
856 * x1' = c[0]*x1 + c[1]*y1 + c[2]
857 * y1' = c[3]*x1 + c[4]*y1 + c[5]
858 * x2' = c[0]*x2 + c[1]*y2 + c[2]
859 * y2' = c[3]*x2 + c[4]*y2 + c[5]
860 * x3' = c[0]*x3 + c[1]*y3 + c[2]
861 * y3' = c[3]*x3 + c[4]*y3 + c[5]
862 *
863 * This can be represented as
864 *
865 * AC = B
866 *
867 * where B and C are column vectors
868 *
869 * B = [ x1' y1' x2' y2' x3' y3' ]
870 * C = [ c[0] c[1] c[2] c[3] c[4] c[5] c[6] ]
871 *
872 * and A is the 6x6 matrix
873 *
874 * x1 y1 1 0 0 0
875 * 0 0 0 x1 y1 1
876 * x2 y2 1 0 0 0
877 * 0 0 0 x2 y2 1
878 * x3 y3 1 0 0 0
879 * 0 0 0 x3 y3 1
880 *
881 * These six equations are solved here for the coefficients C.
882 *
883 * These six coefficients can then be used to find the dest
884 * point x',y') corresponding to any src point (x,y, according
885 * to the equations
886 *
887 * x' = c[0]x + c[1]y + c[2]
888 * y' = c[3]x + c[4]y + c[5]
889 *
890 * that are implemented in affineXformPt.
891 *
892 * !!!!!!!!!!!!!!!!!! Very important !!!!!!!!!!!!!!!!!!!!!!
893 *
894 * When the affine transform is composed from a set of simple
895 * operations such as translation, scaling and rotation,
896 * it is built in a form to convert from the un-transformed src
897 * point to the transformed dest point. However, when an
898 * affine transform is used on images, it is used in an inverted
899 * way: it converts from the transformed dest point to the
900 * un-transformed src point. So, for example, if you transform
901 * a boxa using transform A, to transform an image in the same
902 * way you must use the inverse of A.
903 *
904 * For example, if you transform a boxa with a 3x3 affine matrix
905 * 'mat', the analogous image transformation must use 'matinv':
906 * \code
907 * boxad = boxaAffineTransform(boxas, mat);
908 * affineInvertXform(mat, &matinv);
909 * pixd = pixAffine(pixs, matinv, L_BRING_IN_WHITE);
910 * \endcode
911 * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
912 * </pre>
913 */
914 l_ok
915 getAffineXformCoeffs(PTA *ptas,
916 PTA *ptad,
917 l_float32 **pvc)
918 {
919 l_int32 i;
920 l_float32 x1, y1, x2, y2, x3, y3;
921 l_float32 *b; /* rhs vector of primed coords X'; coeffs returned in *pvc */
922 l_float32 *a[6]; /* 6x6 matrix A */
923
924 if (!ptas)
925 return ERROR_INT("ptas not defined", __func__, 1);
926 if (!ptad)
927 return ERROR_INT("ptad not defined", __func__, 1);
928 if (!pvc)
929 return ERROR_INT("&vc not defined", __func__, 1);
930
931 b = (l_float32 *)LEPT_CALLOC(6, sizeof(l_float32));
932 *pvc = b;
933
934 ptaGetPt(ptas, 0, &x1, &y1);
935 ptaGetPt(ptas, 1, &x2, &y2);
936 ptaGetPt(ptas, 2, &x3, &y3);
937 ptaGetPt(ptad, 0, &b[0], &b[1]);
938 ptaGetPt(ptad, 1, &b[2], &b[3]);
939 ptaGetPt(ptad, 2, &b[4], &b[5]);
940
941 for (i = 0; i < 6; i++)
942 a[i] = (l_float32 *)LEPT_CALLOC(6, sizeof(l_float32));
943 a[0][0] = x1;
944 a[0][1] = y1;
945 a[0][2] = 1.;
946 a[1][3] = x1;
947 a[1][4] = y1;
948 a[1][5] = 1.;
949 a[2][0] = x2;
950 a[2][1] = y2;
951 a[2][2] = 1.;
952 a[3][3] = x2;
953 a[3][4] = y2;
954 a[3][5] = 1.;
955 a[4][0] = x3;
956 a[4][1] = y3;
957 a[4][2] = 1.;
958 a[5][3] = x3;
959 a[5][4] = y3;
960 a[5][5] = 1.;
961
962 gaussjordan(a, b, 6);
963
964 for (i = 0; i < 6; i++)
965 LEPT_FREE(a[i]);
966
967 return 0;
968 }
969
970
971 /*!
972 * \brief affineInvertXform()
973 *
974 * \param[in] vc vector of 6 coefficients
975 * \param[out] pvci inverted transform
976 * \return 0 if OK; 1 on error
977 *
978 * <pre>
979 * Notes:
980 * (1) The 6 affine transform coefficients are the first
981 * two rows of a 3x3 matrix where the last row has
982 * only a 1 in the third column. We invert this
983 * using gaussjordan(), and select the first 2 rows
984 * as the coefficients of the inverse affine transform.
985 * (2) Alternatively, we can find the inverse transform
986 * coefficients by inverting the 2x2 submatrix,
987 * and treating the top 2 coefficients in the 3rd column as
988 * a RHS vector for that 2x2 submatrix. Then the
989 * 6 inverted transform coefficients are composed of
990 * the inverted 2x2 submatrix and the negative of the
991 * transformed RHS vector. Why is this so? We have
992 * Y = AX + R (2 equations in 6 unknowns)
993 * Then
994 * X = A'Y - A'R
995 * Gauss-jordan solves
996 * AF = R
997 * and puts the solution for F, which is A'R,
998 * into the input R vector.
999 *
1000 * </pre>
1001 */
1002 l_ok
1003 affineInvertXform(l_float32 *vc,
1004 l_float32 **pvci)
1005 {
1006 l_int32 i;
1007 l_float32 *vci;
1008 l_float32 *a[3];
1009 l_float32 b[3] = {1.0, 1.0, 1.0}; /* anything; results ignored */
1010
1011 if (!pvci)
1012 return ERROR_INT("&vci not defined", __func__, 1);
1013 *pvci = NULL;
1014 if (!vc)
1015 return ERROR_INT("vc not defined", __func__, 1);
1016
1017 #if 1
1018 for (i = 0; i < 3; i++)
1019 a[i] = (l_float32 *)LEPT_CALLOC(3, sizeof(l_float32));
1020 a[0][0] = vc[0];
1021 a[0][1] = vc[1];
1022 a[0][2] = vc[2];
1023 a[1][0] = vc[3];
1024 a[1][1] = vc[4];
1025 a[1][2] = vc[5];
1026 a[2][2] = 1.0;
1027 gaussjordan(a, b, 3); /* this inverts matrix a */
1028 vci = (l_float32 *)LEPT_CALLOC(6, sizeof(l_float32));
1029 *pvci = vci;
1030 vci[0] = a[0][0];
1031 vci[1] = a[0][1];
1032 vci[2] = a[0][2];
1033 vci[3] = a[1][0];
1034 vci[4] = a[1][1];
1035 vci[5] = a[1][2];
1036 for (i = 0; i < 3; i++)
1037 LEPT_FREE(a[i]);
1038
1039 #else
1040
1041 /* Alternative version, inverting a 2x2 matrix */
1042 { l_float32 *a2[2];
1043 for (i = 0; i < 2; i++)
1044 a2[i] = (l_float32 *)LEPT_CALLOC(2, sizeof(l_float32));
1045 a2[0][0] = vc[0];
1046 a2[0][1] = vc[1];
1047 a2[1][0] = vc[3];
1048 a2[1][1] = vc[4];
1049 b[0] = vc[2];
1050 b[1] = vc[5];
1051 gaussjordan(a2, b, 2); /* this inverts matrix a2 */
1052 vci = (l_float32 *)LEPT_CALLOC(6, sizeof(l_float32));
1053 *pvci = vci;
1054 vci[0] = a2[0][0];
1055 vci[1] = a2[0][1];
1056 vci[2] = -b[0]; /* note sign */
1057 vci[3] = a2[1][0];
1058 vci[4] = a2[1][1];
1059 vci[5] = -b[1]; /* note sign */
1060 for (i = 0; i < 2; i++)
1061 LEPT_FREE(a2[i]);
1062 }
1063 #endif
1064
1065 return 0;
1066 }
1067
1068
1069 /*!
1070 * \brief affineXformSampledPt()
1071 *
1072 * \param[in] vc vector of 6 coefficients
1073 * \param[in] x, y initial point
1074 * \param[out] pxp, pyp transformed point
1075 * \return 0 if OK; 1 on error
1076 *
1077 * <pre>
1078 * Notes:
1079 * (1) This finds the nearest pixel coordinates of the transformed point.
1080 * (2) It does not check ptrs for returned data!
1081 * </pre>
1082 */
1083 l_ok
1084 affineXformSampledPt(l_float32 *vc,
1085 l_int32 x,
1086 l_int32 y,
1087 l_int32 *pxp,
1088 l_int32 *pyp)
1089 {
1090 if (!vc)
1091 return ERROR_INT("vc not defined", __func__, 1);
1092
1093 *pxp = (l_int32)(vc[0] * x + vc[1] * y + vc[2] + 0.5);
1094 *pyp = (l_int32)(vc[3] * x + vc[4] * y + vc[5] + 0.5);
1095 return 0;
1096 }
1097
1098
1099 /*!
1100 * \brief affineXformPt()
1101 *
1102 * \param[in] vc vector of 6 coefficients
1103 * \param[in] x, y initial point
1104 * \param[out] pxp, pyp transformed point
1105 * \return 0 if OK; 1 on error
1106 *
1107 * <pre>
1108 * Notes:
1109 * (1) This computes the floating point location of the transformed point.
1110 * (2) It does not check ptrs for returned data!
1111 * </pre>
1112 */
1113 l_ok
1114 affineXformPt(l_float32 *vc,
1115 l_int32 x,
1116 l_int32 y,
1117 l_float32 *pxp,
1118 l_float32 *pyp)
1119 {
1120 if (!vc)
1121 return ERROR_INT("vc not defined", __func__, 1);
1122
1123 *pxp = vc[0] * x + vc[1] * y + vc[2];
1124 *pyp = vc[3] * x + vc[4] * y + vc[5];
1125 return 0;
1126 }
1127
1128
1129 /*-------------------------------------------------------------*
1130 * Interpolation helper functions *
1131 *-------------------------------------------------------------*/
1132 /*!
1133 * \brief linearInterpolatePixelColor()
1134 *
1135 * \param[in] datas ptr to beginning of image data
1136 * \param[in] wpls 32-bit word/line for this data array
1137 * \param[in] w, h of image
1138 * \param[in] x, y floating pt location for evaluation
1139 * \param[in] colorval color brought in from the outside when the
1140 * input x,y location is outside the image;
1141 * in 0xrrggbb00 format)
1142 * \param[out] pval interpolated color value
1143 * \return 0 if OK, 1 on error
1144 *
1145 * <pre>
1146 * Notes:
1147 * (1) This is a standard linear interpolation function. It is
1148 * equivalent to area weighting on each component, and
1149 * avoids "jaggies" when rendering sharp edges.
1150 * </pre>
1151 */
1152 l_ok
1153 linearInterpolatePixelColor(l_uint32 *datas,
1154 l_int32 wpls,
1155 l_int32 w,
1156 l_int32 h,
1157 l_float32 x,
1158 l_float32 y,
1159 l_uint32 colorval,
1160 l_uint32 *pval)
1161 {
1162 l_int32 valid, xpm, ypm, xp, xp2, yp, xf, yf;
1163 l_int32 rval, gval, bval;
1164 l_uint32 word00, word01, word10, word11;
1165 l_uint32 *lines;
1166
1167 if (!pval)
1168 return ERROR_INT("&val not defined", __func__, 1);
1169 *pval = colorval;
1170 if (!datas)
1171 return ERROR_INT("datas not defined", __func__, 1);
1172
1173 /* Skip if x or y are invalid. (x,y) must be in the source image.
1174 * Failure to detect an invalid point will cause a mem address fault.
1175 * Occasionally, x or y will be a nan, and relational checks always
1176 * fail for nans. Therefore we check if the point is inside the pix */
1177 valid = (x >= 0.0 && y >= 0.0 && x < w && y < h);
1178 if (!valid) return 0;
1179
1180 xpm = (l_int32)(16.0 * x);
1181 ypm = (l_int32)(16.0 * y);
1182 xp = xpm >> 4;
1183 xp2 = xp + 1 < w ? xp + 1 : xp;
1184 yp = ypm >> 4;
1185 if (yp + 1 >= h) wpls = 0;
1186 xf = xpm & 0x0f;
1187 yf = ypm & 0x0f;
1188
1189 #if DEBUG
1190 if (xf < 0 || yf < 0)
1191 lept_stderr("xp = %d, yp = %d, xf = %d, yf = %d\n", xp, yp, xf, yf);
1192 #endif /* DEBUG */
1193
1194 /* Do area weighting (eqiv. to linear interpolation) */
1195 lines = datas + yp * wpls;
1196 word00 = *(lines + xp);
1197 word10 = *(lines + xp2);
1198 word01 = *(lines + wpls + xp);
1199 word11 = *(lines + wpls + xp2);
1200 rval = ((16 - xf) * (16 - yf) * ((word00 >> L_RED_SHIFT) & 0xff) +
1201 xf * (16 - yf) * ((word10 >> L_RED_SHIFT) & 0xff) +
1202 (16 - xf) * yf * ((word01 >> L_RED_SHIFT) & 0xff) +
1203 xf * yf * ((word11 >> L_RED_SHIFT) & 0xff)) / 256;
1204 gval = ((16 - xf) * (16 - yf) * ((word00 >> L_GREEN_SHIFT) & 0xff) +
1205 xf * (16 - yf) * ((word10 >> L_GREEN_SHIFT) & 0xff) +
1206 (16 - xf) * yf * ((word01 >> L_GREEN_SHIFT) & 0xff) +
1207 xf * yf * ((word11 >> L_GREEN_SHIFT) & 0xff)) / 256;
1208 bval = ((16 - xf) * (16 - yf) * ((word00 >> L_BLUE_SHIFT) & 0xff) +
1209 xf * (16 - yf) * ((word10 >> L_BLUE_SHIFT) & 0xff) +
1210 (16 - xf) * yf * ((word01 >> L_BLUE_SHIFT) & 0xff) +
1211 xf * yf * ((word11 >> L_BLUE_SHIFT) & 0xff)) / 256;
1212 composeRGBPixel(rval, gval, bval, pval);
1213 return 0;
1214 }
1215
1216
1217 /*!
1218 * \brief linearInterpolatePixelGray()
1219 *
1220 * \param[in] datas ptr to beginning of image data
1221 * \param[in] wpls 32-bit word/line for this data array
1222 * \param[in] w, h of image
1223 * \param[in] x, y floating pt location for evaluation
1224 * \param[in] grayval color brought in from the outside when the
1225 * input x,y location is outside the image
1226 * \param[out] pval interpolated gray value
1227 * \return 0 if OK, 1 on error
1228 *
1229 * <pre>
1230 * Notes:
1231 * (1) This is a standard linear interpolation function. It is
1232 * equivalent to area weighting on each component, and
1233 * avoids "jaggies" when rendering sharp edges.
1234 * </pre>
1235 */
1236 l_ok
1237 linearInterpolatePixelGray(l_uint32 *datas,
1238 l_int32 wpls,
1239 l_int32 w,
1240 l_int32 h,
1241 l_float32 x,
1242 l_float32 y,
1243 l_int32 grayval,
1244 l_int32 *pval)
1245 {
1246 l_int32 valid, xpm, ypm, xp, xp2, yp, xf, yf, v00, v10, v01, v11;
1247 l_uint32 *lines;
1248
1249 if (!pval)
1250 return ERROR_INT("&val not defined", __func__, 1);
1251 *pval = grayval;
1252 if (!datas)
1253 return ERROR_INT("datas not defined", __func__, 1);
1254
1255 /* Skip if x or y is invalid. (x,y) must be in the source image.
1256 * Failure to detect an invalid point will cause a mem address fault.
1257 * Occasionally, x or y will be a nan, and relational checks always
1258 * fail for nans. Therefore we check if the point is inside the pix */
1259 valid = (x >= 0.0 && y >= 0.0 && x < w && y < h);
1260 if (!valid) return 0;
1261
1262 xpm = (l_int32)(16.0 * x);
1263 ypm = (l_int32)(16.0 * y);
1264 xp = xpm >> 4;
1265 xp2 = xp + 1 < w ? xp + 1 : xp;
1266 yp = ypm >> 4;
1267 if (yp + 1 >= h) wpls = 0;
1268 xf = xpm & 0x0f;
1269 yf = ypm & 0x0f;
1270
1271 #if DEBUG
1272 if (xf < 0 || yf < 0)
1273 lept_stderr("xp = %d, yp = %d, xf = %d, yf = %d\n", xp, yp, xf, yf);
1274 #endif /* DEBUG */
1275
1276 /* Interpolate by area weighting. */
1277 lines = datas + yp * wpls;
1278 v00 = (16 - xf) * (16 - yf) * GET_DATA_BYTE(lines, xp);
1279 v10 = xf * (16 - yf) * GET_DATA_BYTE(lines, xp2);
1280 v01 = (16 - xf) * yf * GET_DATA_BYTE(lines + wpls, xp);
1281 v11 = xf * yf * GET_DATA_BYTE(lines + wpls, xp2);
1282 *pval = (v00 + v01 + v10 + v11) / 256;
1283 return 0;
1284 }
1285
1286
1287
1288 /*-------------------------------------------------------------*
1289 * Gauss-jordan linear equation solver *
1290 *-------------------------------------------------------------*/
1291 #define SWAP(a,b) {temp = (a); (a) = (b); (b) = temp;}
1292
1293 /*!
1294 * \brief gaussjordan()
1295 *
1296 * \param[in] a n x n matrix
1297 * \param[in] b n x 1 right-hand side column vector
1298 * \param[in] n dimension
1299 * \return 0 if ok, 1 on error
1300 *
1301 * <pre>
1302 * Notes:
1303 * (1) There are two side-effects:
1304 * * The matrix a is transformed to its inverse A
1305 * * The rhs vector b is transformed to the solution x
1306 * of the linear equation ax = b
1307 * (2) The inverse A can then be used to solve the same equation with
1308 * different rhs vectors c by multiplication: x = Ac
1309 * (3) Adapted from "Numerical Recipes in C, Second Edition", 1992,
1310 * pp. 36-41 (gauss-jordan elimination)
1311 * </pre>
1312 */
1313 l_int32
1314 gaussjordan(l_float32 **a,
1315 l_float32 *b,
1316 l_int32 n)
1317 {
1318 l_int32 i, icol, irow, j, k, col, row, success;
1319 l_int32 *indexc, *indexr, *ipiv;
1320 l_float32 maxval, val, pivinv, temp;
1321
1322 if (!a)
1323 return ERROR_INT("a not defined", __func__, 1);
1324 if (!b)
1325 return ERROR_INT("b not defined", __func__, 1);
1326
1327 success = TRUE;
1328 indexc = (l_int32 *)LEPT_CALLOC(n, sizeof(l_int32));
1329 indexr = (l_int32 *)LEPT_CALLOC(n, sizeof(l_int32));
1330 ipiv = (l_int32 *)LEPT_CALLOC(n, sizeof(l_int32));
1331 if (!indexc || !indexr || !ipiv) {
1332 L_ERROR("array not made\n", __func__);
1333 success = FALSE;
1334 goto cleanup_arrays;
1335 }
1336
1337 icol = irow = 0; /* silence static checker */
1338 for (i = 0; i < n; i++) {
1339 maxval = 0.0;
1340 for (j = 0; j < n; j++) {
1341 if (ipiv[j] != 1) {
1342 for (k = 0; k < n; k++) {
1343 if (ipiv[k] == 0) {
1344 if (fabs(a[j][k]) >= maxval) {
1345 maxval = fabs(a[j][k]);
1346 irow = j;
1347 icol = k;
1348 }
1349 } else if (ipiv[k] > 1) {
1350 L_ERROR("singular matrix\n", __func__);
1351 success = FALSE;
1352 goto cleanup_arrays;
1353 }
1354 }
1355 }
1356 }
1357 ++(ipiv[icol]);
1358
1359 if (irow != icol) {
1360 for (col = 0; col < n; col++)
1361 SWAP(a[irow][col], a[icol][col]);
1362 SWAP(b[irow], b[icol]);
1363 }
1364
1365 indexr[i] = irow;
1366 indexc[i] = icol;
1367 if (a[icol][icol] == 0.0) {
1368 L_ERROR("singular matrix\n", __func__);
1369 success = FALSE;
1370 goto cleanup_arrays;
1371 }
1372 pivinv = 1.0 / a[icol][icol];
1373 a[icol][icol] = 1.0;
1374 for (col = 0; col < n; col++)
1375 a[icol][col] *= pivinv;
1376 b[icol] *= pivinv;
1377
1378 for (row = 0; row < n; row++) {
1379 if (row != icol) {
1380 val = a[row][icol];
1381 a[row][icol] = 0.0;
1382 for (col = 0; col < n; col++)
1383 a[row][col] -= a[icol][col] * val;
1384 b[row] -= b[icol] * val;
1385 }
1386 }
1387 }
1388
1389 for (col = n - 1; col >= 0; col--) {
1390 if (indexr[col] != indexc[col]) {
1391 for (k = 0; k < n; k++)
1392 SWAP(a[k][indexr[col]], a[k][indexc[col]]);
1393 }
1394 }
1395
1396 cleanup_arrays:
1397 LEPT_FREE(indexr);
1398 LEPT_FREE(indexc);
1399 LEPT_FREE(ipiv);
1400 return (success) ? 0 : 1;
1401 }
1402
1403
1404 /*-------------------------------------------------------------*
1405 * Sequential affine image transformation *
1406 *-------------------------------------------------------------*/
1407 /*!
1408 * \brief pixAffineSequential()
1409 *
1410 * \param[in] pixs
1411 * \param[in] ptad 3 pts of final coordinate space
1412 * \param[in] ptas 3 pts of initial coordinate space
1413 * \param[in] bw pixels of additional border width during computation
1414 * \param[in] bh pixels of additional border height during computation
1415 * \return pixd, or NULL on error
1416 *
1417 * <pre>
1418 * Notes:
1419 * (1) The 3 pts must not be collinear.
1420 * (2) The 3 pts must be given in this order:
1421 * ~ origin
1422 * ~ a location along the x-axis
1423 * ~ a location along the y-axis.
1424 * (3) You must guess how much border must be added so that no
1425 * pixels are lost in the transformations from src to
1426 * dest coordinate space. (This can be calculated but it
1427 * is a lot of work!) For coordinate spaces that are nearly
1428 * at right angles, on a 300 ppi scanned page, the addition
1429 * of 1000 pixels on each side is usually sufficient.
1430 * (4) This is here for pedagogical reasons. It is about 3x faster
1431 * on 1 bpp images than pixAffineSampled(), but the results
1432 * on text are much inferior.
1433 * </pre>
1434 */
1435 PIX *
1436 pixAffineSequential(PIX *pixs,
1437 PTA *ptad,
1438 PTA *ptas,
1439 l_int32 bw,
1440 l_int32 bh)
1441 {
1442 l_int32 x1, y1, x2, y2, x3, y3; /* ptas */
1443 l_int32 x1p, y1p, x2p, y2p, x3p, y3p; /* ptad */
1444 l_int32 x1sc, y1sc; /* scaled origin */
1445 l_float32 x2s, x2sp, scalex, scaley;
1446 l_float32 th3, th3p, ph2, ph2p;
1447 #if DEBUG
1448 l_float32 rad2deg;
1449 #endif /* DEBUG */
1450 PIX *pix1, *pix2, *pixd;
1451
1452 if (!pixs)
1453 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
1454 if (!ptas)
1455 return (PIX *)ERROR_PTR("ptas not defined", __func__, NULL);
1456 if (!ptad)
1457 return (PIX *)ERROR_PTR("ptad not defined", __func__, NULL);
1458
1459 if (ptaGetCount(ptas) != 3)
1460 return (PIX *)ERROR_PTR("ptas count not 3", __func__, NULL);
1461 if (ptaGetCount(ptad) != 3)
1462 return (PIX *)ERROR_PTR("ptad count not 3", __func__, NULL);
1463 ptaGetIPt(ptas, 0, &x1, &y1);
1464 ptaGetIPt(ptas, 1, &x2, &y2);
1465 ptaGetIPt(ptas, 2, &x3, &y3);
1466 ptaGetIPt(ptad, 0, &x1p, &y1p);
1467 ptaGetIPt(ptad, 1, &x2p, &y2p);
1468 ptaGetIPt(ptad, 2, &x3p, &y3p);
1469
1470 pix1 = pix2 = pixd = NULL;
1471
1472 if (y1 == y3)
1473 return (PIX *)ERROR_PTR("y1 == y3!", __func__, NULL);
1474 if (y1p == y3p)
1475 return (PIX *)ERROR_PTR("y1p == y3p!", __func__, NULL);
1476
1477 if (bw != 0 || bh != 0) {
1478 /* resize all points and add border to pixs */
1479 x1 = x1 + bw;
1480 y1 = y1 + bh;
1481 x2 = x2 + bw;
1482 y2 = y2 + bh;
1483 x3 = x3 + bw;
1484 y3 = y3 + bh;
1485 x1p = x1p + bw;
1486 y1p = y1p + bh;
1487 x2p = x2p + bw;
1488 y2p = y2p + bh;
1489 x3p = x3p + bw;
1490 y3p = y3p + bh;
1491
1492 if ((pix1 = pixAddBorderGeneral(pixs, bw, bw, bh, bh, 0)) == NULL)
1493 return (PIX *)ERROR_PTR("pix1 not made", __func__, NULL);
1494 } else {
1495 pix1 = pixCopy(NULL, pixs);
1496 }
1497
1498 /*-------------------------------------------------------------*
1499 The horizontal shear is done to move the 3rd point to the
1500 y axis. This moves the 2nd point either towards or away
1501 from the y axis, depending on whether it is above or below
1502 the x axis. That motion must be computed so that we know
1503 the angle of vertical shear to use to get the 2nd point
1504 on the x axis. We must also know the x coordinate of the
1505 2nd point in order to compute how much scaling is required
1506 to match points on the axis.
1507 *-------------------------------------------------------------*/
1508
1509 /* Shear angles required to put src points on x and y axes */
1510 th3 = atan2((l_float64)(x1 - x3), (l_float64)(y1 - y3));
1511 x2s = (l_float32)(x2 - ((l_float32)(y1 - y2) * (x3 - x1)) / (y1 - y3));
1512 if (x2s == (l_float32)x1) {
1513 L_ERROR("x2s == x1!\n", __func__);
1514 goto cleanup_pix;
1515 }
1516 ph2 = atan2((l_float64)(y1 - y2), (l_float64)(x2s - x1));
1517
1518 /* Shear angles required to put dest points on x and y axes.
1519 * Use the negative of these values to instead move the
1520 * src points from the axes to the actual dest position.
1521 * These values are also needed to scale the image. */
1522 th3p = atan2((l_float64)(x1p - x3p), (l_float64)(y1p - y3p));
1523 x2sp = (l_float32)(x2p -
1524 ((l_float32)(y1p - y2p) * (x3p - x1p)) / (y1p - y3p));
1525 if (x2sp == (l_float32)x1p) {
1526 L_ERROR("x2sp == x1p!\n", __func__);
1527 goto cleanup_pix;
1528 }
1529 ph2p = atan2((l_float64)(y1p - y2p), (l_float64)(x2sp - x1p));
1530
1531 /* Shear image to first put src point 3 on the y axis,
1532 * and then to put src point 2 on the x axis */
1533 pixHShearIP(pix1, y1, th3, L_BRING_IN_WHITE);
1534 pixVShearIP(pix1, x1, ph2, L_BRING_IN_WHITE);
1535
1536 /* Scale image to match dest scale. The dest scale
1537 * is calculated above from the angles th3p and ph2p
1538 * that would be required to move the dest points to
1539 * the x and y axes. */
1540 scalex = (l_float32)(x2sp - x1p) / (x2s - x1);
1541 scaley = (l_float32)(y3p - y1p) / (y3 - y1);
1542 if ((pix2 = pixScale(pix1, scalex, scaley)) == NULL) {
1543 L_ERROR("pix2 not made\n", __func__);
1544 goto cleanup_pix;
1545 }
1546
1547 #if DEBUG
1548 rad2deg = 180. / 3.1415926535;
1549 lept_stderr("th3 = %5.1f deg, ph2 = %5.1f deg\n",
1550 rad2deg * th3, rad2deg * ph2);
1551 lept_stderr("th3' = %5.1f deg, ph2' = %5.1f deg\n",
1552 rad2deg * th3p, rad2deg * ph2p);
1553 lept_stderr("scalex = %6.3f, scaley = %6.3f\n", scalex, scaley);
1554 #endif /* DEBUG */
1555
1556 /*-------------------------------------------------------------*
1557 Scaling moves the 1st src point, which is the origin.
1558 It must now be moved again to coincide with the origin
1559 (1st point) of the dest. After this is done, the 2nd
1560 and 3rd points must be sheared back to the original
1561 positions of the 2nd and 3rd dest points. We use the
1562 negative of the angles that were previously computed
1563 for shearing those points in the dest image to x and y
1564 axes, and take the shears in reverse order as well.
1565 *-------------------------------------------------------------*/
1566 /* Shift image to match dest origin. */
1567 x1sc = (l_int32)(scalex * x1 + 0.5); /* x comp of origin after scaling */
1568 y1sc = (l_int32)(scaley * y1 + 0.5); /* y comp of origin after scaling */
1569 pixRasteropIP(pix2, x1p - x1sc, y1p - y1sc, L_BRING_IN_WHITE);
1570
1571 /* Shear image to take points 2 and 3 off the axis and
1572 * put them in the original dest position */
1573 pixVShearIP(pix2, x1p, -ph2p, L_BRING_IN_WHITE);
1574 pixHShearIP(pix2, y1p, -th3p, L_BRING_IN_WHITE);
1575
1576 if (bw != 0 || bh != 0) {
1577 if ((pixd = pixRemoveBorderGeneral(pix2, bw, bw, bh, bh)) == NULL)
1578 L_ERROR("pixd not made\n", __func__);
1579 } else {
1580 pixd = pixClone(pix2);
1581 }
1582
1583 cleanup_pix:
1584 pixDestroy(&pix1);
1585 pixDestroy(&pix2);
1586 return pixd;
1587 }