/*====================================================================*
- Copyright (C) 2001 Leptonica. All rights reserved.
- This software is distributed in the hope that it will be
- useful, but with NO WARRANTY OF ANY KIND.
- No author or distributor accepts responsibility to anyone for the
- consequences of using this software, or for whether it serves any
- particular purpose or works at all, unless he or she says so in
- writing. Everyone is granted permission to copy, modify and
- redistribute this source code, for commercial or non-commercial
- purposes, with the following restrictions: (1) the origin of this
- source code must not be misrepresented; (2) modified versions must
- be plainly marked as such; and (3) this notice may not be removed
- or altered from any source or modified source distribution.
*====================================================================*/
/*
* affinecompose.c
*
* Composable coordinate transforms
* l_float32 *createMatrix2dTranslate()
* l_float32 *createMatrixScale()
* l_float32 *createMatrixRotate()
*
* Special coordinate transforms on pta
* PTA *ptaTranslate()
* PTA *ptaScale()
* PTA *ptaRotate()
*
* Special coordinate transforms on boxa
* BOXA *boxaTranslate()
* BOXA *boxaScale()
* BOXA *boxaRotate()
*
* General coordinate transform on pta and boxa
* PTA *ptaAffineTransform()
* BOXA *boxaAffineTransform()
*
* Matrix operations
* l_int32 l_productMatVec()
* l_int32 l_productMat2()
* l_int32 l_productMat3()
* l_int32 l_productMat4()
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "allheaders.h"
/*-------------------------------------------------------------*
* Composable coordinate transforms *
*-------------------------------------------------------------*/
/*!
* createMatrix2dTranslate()
*
* Input: transx (x component of translation wrt. the origin)
* transy (y component of translation wrt. the origin)
* Return: 3x3 transform matrix, or null on error
*
* Notes;
* (1) The translation is equivalent to:
* v' = Av
* where v and v' are 1x3 column vectors in the form
* v = [x, y, 1]^ (^ denotes transpose)
* and the affine tranlation matrix is
* A = [ 1 0 tx
* 0 1 ty
* 0 0 1 ]
*
* (2) We consider translation as with respect to a fixed origin.
* In a clipping operation, the origin moves and the points
* are fixed, and you use (-tx, -ty) where (tx, ty) is the
* translation vector of the origin.
*/
l_float32 *
createMatrix2dTranslate(l_float32 transx,
l_float32 transy)
{
l_float32 *mat;
PROCNAME("createMatrix2dTranslate");
if ((mat = (l_float32 *)CALLOC(9, sizeof(l_float32))) == NULL)
return (l_float32 *)ERROR_PTR("mat not made", procName, NULL);
mat[0] = mat[4] = mat[8] = 1;
mat[2] = transx;
mat[5] = transy;
return mat;
}
/*!
* createMatrix2dScale()
*
* Input: scalex (horizontal scale factor)
* scaley (vertical scale factor)
* Return: 3x3 transform matrix, or null on error
*
* Notes;
* (1) The scaling is equivalent to:
* v' = Av
* where v and v' are 1x3 column vectors in the form
* v = [x, y, 1]^ (^ denotes transpose)
* and the affine scaling matrix is
* A = [ sx 0 0
* 0 sy 0
* 0 0 1 ]
*
* (2) We consider scaling as with respect to a fixed origin.
* In other words, the origin is the only point that doesn't
* move in the scaling transform.
*/
l_float32 *
createMatrix2dScale(l_float32 scalex,
l_float32 scaley)
{
l_float32 *mat;
PROCNAME("createMatrix2dScale");
if ((mat = (l_float32 *)CALLOC(9, sizeof(l_float32))) == NULL)
return (l_float32 *)ERROR_PTR("mat not made", procName, NULL);
mat[0] = scalex;
mat[4] = scaley;
mat[8] = 1;
return mat;
}
/*!
* createMatrix2dRotate()
*
* Input: xc, yc (location of center of rotation)
* angle (rotation in radians; clockwise is positive)
* Return: 3x3 transform matrix, or null on error
*
* Notes;
* (1) The rotation is equivalent to:
* v' = Av
* where v and v' are 1x3 column vectors in the form
* v = [x, y, 1]^ (^ denotes transpose)
* and the affine rotation matrix is
* A = [ cosa -sina xc*(1-cosa) + yc*sina
* sina cosa yc*(1-cosa) - xc*sina
* 0 0 1 ]
*
* If the rotation is about the origin, (xc, yc) = (0, 0) and
* this simplifies to
* A = [ cosa -sina 0
* sina cosa 0
* 0 0 1 ]
*
* These relations follow from the following equations, which
* you can convince yourself are correct as follows. Draw a
* circle centered on (xc,yc) and passing through (x,y), with
* (x',y') on the arc at an angle 'a' clockwise from (x,y).
* [ Hint: cos(a + b) = cosa * cosb - sina * sinb
* sin(a + b) = sina * cosb + cosa * sinb ]
*
* x' - xc = (x - xc) * cosa - (y - yc) * sina
* y' - yc = (x - xc) * sina + (y - yc) * cosa
*/
l_float32 *
createMatrix2dRotate(l_float32 xc,
l_float32 yc,
l_float32 angle)
{
l_float32 sina, cosa;
l_float32 *mat;
PROCNAME("createMatrix2dRotate");
if ((mat = (l_float32 *)CALLOC(9, sizeof(l_float32))) == NULL)
return (l_float32 *)ERROR_PTR("mat not made", procName, NULL);
sina = sin(angle);
cosa = cos(angle);
mat[0] = mat[4] = cosa;
mat[1] = -sina;
mat[2] = xc * (1.0 - cosa) + yc * sina;
mat[3] = sina;
mat[5] = yc * (1.0 - cosa) - xc * sina;
mat[8] = 1;
return mat;
}
/*-------------------------------------------------------------*
* Special coordinate transforms on pta *
*-------------------------------------------------------------*/
/*!
* ptaTranslate()
*
* Input: ptas (for initial points)
* transx (x component of translation wrt. the origin)
* transy (y component of translation wrt. the origin)
* Return: ptad (translated points), or null on error
*
* Notes;
* (1) See createMatrix2dTranslate() for details of transform.
*/
PTA *
ptaTranslate(PTA *ptas,
l_float32 transx,
l_float32 transy)
{
l_int32 i, npts;
l_float32 x, y;
PTA *ptad;
PROCNAME("ptaTranslate");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
npts = ptaGetCount(ptas);
if ((ptad = ptaCreate(npts)) == NULL)
return (PTA *)ERROR_PTR("ptad not made", procName, NULL);
for (i = 0; i < npts; i++) {
ptaGetPt(ptas, i, &x, &y);
ptaAddPt(ptad, x + transx, y + transy);
}
return ptad;
}
/*!
* ptaScale()
*
* Input: ptas (for initial points)
* scalex (horizontal scale factor)
* scaley (vertical scale factor)
* Return: 0 if OK; 1 on error
*
* Notes;
* (1) See createMatrix2dScale() for details of transform.
*/
PTA *
ptaScale(PTA *ptas,
l_float32 scalex,
l_float32 scaley)
{
l_int32 i, npts;
l_float32 x, y;
PTA *ptad;
PROCNAME("ptaScale");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
npts = ptaGetCount(ptas);
if ((ptad = ptaCreate(npts)) == NULL)
return (PTA *)ERROR_PTR("ptad not made", procName, NULL);
for (i = 0; i < npts; i++) {
ptaGetPt(ptas, i, &x, &y);
ptaAddPt(ptad, scalex * x, scaley * y);
}
return ptad;
}
/*!
* ptaRotate()
*
* Input: ptas (for initial points)
* (xc, yc) (location of center of rotation)
* angle (rotation in radians; clockwise is positive)
* (&ptad) (<return> new locations)
* Return: 0 if OK; 1 on error
*
* Notes;
* (1) See createMatrix2dScale() for details of transform.
*/
PTA *
ptaRotate(PTA *ptas,
l_float32 xc,
l_float32 yc,
l_float32 angle)
{
l_int32 i, npts;
l_float32 x, y, xp, yp, sina, cosa;
PTA *ptad;
PROCNAME("ptaRotate");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
npts = ptaGetCount(ptas);
if ((ptad = ptaCreate(npts)) == NULL)
return (PTA *)ERROR_PTR("ptad not made", procName, NULL);
sina = sin(angle);
cosa = cos(angle);
for (i = 0; i < npts; i++) {
ptaGetPt(ptas, i, &x, &y);
xp = xc + (x - xc) * cosa - (y - yc) * sina;
yp = yc + (x - xc) * sina + (y - yc) * cosa;
ptaAddPt(ptad, xp, yp);
}
return ptad;
}
/*-------------------------------------------------------------*
* Special coordinate transforms on boxa *
*-------------------------------------------------------------*/
/*!
* boxaTranslate()
*
* Input: boxas
* transx (x component of translation wrt. the origin)
* transy (y component of translation wrt. the origin)
* Return: boxad (translated boxas), or null on error
*
* Notes;
* (1) See createMatrix2dTranslate() for details of transform.
*/
BOXA *
boxaTranslate(BOXA *boxas,
l_float32 transx,
l_float32 transy)
{
PTA *ptas, *ptad;
BOXA *boxad;
PROCNAME("boxaTranslate");
if (!boxas)
return (BOXA *)ERROR_PTR("boxas not defined", procName, NULL);
ptas = boxaConvertToPta(boxas, 4);
ptad = ptaTranslate(ptas, transx, transy);
boxad = ptaConvertToBoxa(ptad, 4);
ptaDestroy(&ptas);
ptaDestroy(&ptad);
return boxad;
}
/*!
* boxaScale()
*
* Input: boxas
* scalex (horizontal scale factor)
* scaley (vertical scale factor)
* Return: boxad (scaled boxas), or null on error
*
* Notes;
* (1) See createMatrix2dScale() for details of transform.
*/
BOXA *
boxaScale(BOXA *boxas,
l_float32 scalex,
l_float32 scaley)
{
PTA *ptas, *ptad;
BOXA *boxad;
PROCNAME("boxaScale");
if (!boxas)
return (BOXA *)ERROR_PTR("boxas not defined", procName, NULL);
ptas = boxaConvertToPta(boxas, 4);
ptad = ptaScale(ptas, scalex, scaley);
boxad = ptaConvertToBoxa(ptad, 4);
ptaDestroy(&ptas);
ptaDestroy(&ptad);
return boxad;
}
/*!
* boxaRotate()
*
* Input: boxas
* (xc, yc) (location of center of rotation)
* angle (rotation in radians; clockwise is positive)
* Return: boxad (scaled boxas), or null on error
*
* Notes;
* (1) See createMatrix2dRotate() for details of transform.
*/
BOXA *
boxaRotate(BOXA *boxas,
l_float32 xc,
l_float32 yc,
l_float32 angle)
{
PTA *ptas, *ptad;
BOXA *boxad;
PROCNAME("boxaRotate");
if (!boxas)
return (BOXA *)ERROR_PTR("boxas not defined", procName, NULL);
ptas = boxaConvertToPta(boxas, 4);
ptad = ptaRotate(ptas, xc, yc, angle);
boxad = ptaConvertToBoxa(ptad, 4);
ptaDestroy(&ptas);
ptaDestroy(&ptad);
return boxad;
}
/*-------------------------------------------------------------*
* General affine coordinate transform *
*-------------------------------------------------------------*/
/*!
* ptaAffineTransform()
*
* Input: ptas (for initial points)
* mat (3x3 transform matrix; canonical form)
* Return: ptad (transformed points), or null on error
*/
PTA *
ptaAffineTransform(PTA *ptas,
l_float32 *mat)
{
l_int32 i, npts;
l_float32 vecs[3], vecd[3];
PTA *ptad;
PROCNAME("ptaAffineTransform");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
if (!mat)
return (PTA *)ERROR_PTR("transform not defined", procName, NULL);
vecs[2] = 1;
npts = ptaGetCount(ptas);
if ((ptad = ptaCreate(npts)) == NULL)
return (PTA *)ERROR_PTR("ptad not made", procName, NULL);
for (i = 0; i < npts; i++) {
ptaGetPt(ptas, i, &vecs[0], &vecs[1]);
l_productMatVec(mat, vecs, vecd, 3);
ptaAddPt(ptad, vecd[0], vecd[1]);
}
return ptad;
}
/*!
* boxaAffineTransform()
*
* Input: boxas
* mat (3x3 transform matrix; canonical form)
* Return: boxad (transformed boxas), or null on error
*/
BOXA *
boxaAffineTransform(BOXA *boxas,
l_float32 *mat)
{
PTA *ptas, *ptad;
BOXA *boxad;
PROCNAME("boxaAffineTransform");
if (!boxas)
return (BOXA *)ERROR_PTR("boxas not defined", procName, NULL);
if (!mat)
return (BOXA *)ERROR_PTR("transform not defined", procName, NULL);
ptas = boxaConvertToPta(boxas, 4);
ptad = ptaAffineTransform(ptas, mat);
boxad = ptaConvertToBoxa(ptad, 4);
ptaDestroy(&ptas);
ptaDestroy(&ptad);
return boxad;
}
/*-------------------------------------------------------------*
* Matrix operations *
*-------------------------------------------------------------*/
/*!
* l_productMatVec()
*
* Input: mat (square matrix, as a 1-dimensional @size^2 array)
* vecs (input column vector of length @size)
* vecd (result column vector)
* size (matrix is @size x @size; vectors are length @size)
* Return: 0 if OK, 1 on error
*/
l_int32
l_productMatVec(l_float32 *mat,
l_float32 *vecs,
l_float32 *vecd,
l_int32 size)
{
l_int32 i, j;
PROCNAME("l_productMatVec");
if (!mat)
return ERROR_INT("matrix not defined", procName, 1);
if (!vecs)
return ERROR_INT("input vector not defined", procName, 1);
if (!vecd)
return ERROR_INT("result vector not defined", procName, 1);
for (i = 0; i < size; i++) {
vecd[i] = 0;
for (j = 0; j < size; j++) {
vecd[i] += mat[size * i + j] * vecs[j];
}
}
return 0;
}
/*!
* l_productMat2()
*
* Input: mat1 (square matrix, as a 1-dimensional size^2 array)
* mat2 (square matrix, as a 1-dimensional size^2 array)
* matd (square matrix; product stored here)
* size (of matrices)
* Return: 0 if OK, 1 on error
*/
l_int32
l_productMat2(l_float32 *mat1,
l_float32 *mat2,
l_float32 *matd,
l_int32 size)
{
l_int32 i, j, k, index;
PROCNAME("l_productMat2");
if (!mat1)
return ERROR_INT("matrix 1 not defined", procName, 1);
if (!mat2)
return ERROR_INT("matrix 2 not defined", procName, 1);
if (!matd)
return ERROR_INT("result matrix not defined", procName, 1);
for (i = 0; i < size; i++) {
for (j = 0; j < size; j++) {
index = size * i + j;
matd[index] = 0;
for (k = 0; k < size; k++)
matd[index] += mat1[size * i + k] * mat2[size * k + j];
}
}
return 0;
}
/*!
* l_productMat3()
*
* Input: mat1 (square matrix, as a 1-dimensional size^2 array)
* mat2 (square matrix, as a 1-dimensional size^2 array)
* mat3 (square matrix, as a 1-dimensional size^2 array)
* matd (square matrix; product stored here)
* size (of matrices)
* Return: 0 if OK, 1 on error
*/
l_int32
l_productMat3(l_float32 *mat1,
l_float32 *mat2,
l_float32 *mat3,
l_float32 *matd,
l_int32 size)
{
l_float32 *matt;
PROCNAME("l_productMat3");
if (!mat1)
return ERROR_INT("matrix 1 not defined", procName, 1);
if (!mat2)
return ERROR_INT("matrix 2 not defined", procName, 1);
if (!mat3)
return ERROR_INT("matrix 3 not defined", procName, 1);
if (!matd)
return ERROR_INT("result matrix not defined", procName, 1);
if ((matt = (l_float32 *)CALLOC(size * size, sizeof(l_float32))) == NULL)
return ERROR_INT("matt not made", procName, 1);
l_productMat2(mat1, mat2, matt, size);
l_productMat2(matt, mat3, matd, size);
FREE(matt);
return 0;
}
/*!
* l_productMat4()
*
* Input: mat1 (square matrix, as a 1-dimensional size^2 array)
* mat2 (square matrix, as a 1-dimensional size^2 array)
* mat3 (square matrix, as a 1-dimensional size^2 array)
* mat4 (square matrix, as a 1-dimensional size^2 array)
* matd (square matrix; product stored here)
* size (of matrices)
* Return: 0 if OK, 1 on error
*/
l_int32
l_productMat4(l_float32 *mat1,
l_float32 *mat2,
l_float32 *mat3,
l_float32 *mat4,
l_float32 *matd,
l_int32 size)
{
l_float32 *matt;
PROCNAME("l_productMat4");
if (!mat1)
return ERROR_INT("matrix 1 not defined", procName, 1);
if (!mat2)
return ERROR_INT("matrix 2 not defined", procName, 1);
if (!mat3)
return ERROR_INT("matrix 3 not defined", procName, 1);
if (!matd)
return ERROR_INT("result matrix not defined", procName, 1);
if ((matt = (l_float32 *)CALLOC(size * size, sizeof(l_float32))) == NULL)
return ERROR_INT("matt not made", procName, 1);
l_productMat3(mat1, mat2, mat3, matt, size);
l_productMat2(matt, mat4, matd, size);
FREE(matt);
return 0;
}