sinosoidal

/*
 * pmap_quad_rect: find mapping between quadrilateral and rectangle.
 * The correspondence is:
 *
 *      quad[0] --> (u0,v0)
 *      quad[1] --> (u1,v0)
 *      quad[2] --> (u1,v1)
 *      quad[3] --> (u0,v1)
 *
 * This method of computing the adjoint numerically is cheaper than
 * computing it symbolically.
 */

int pmap_quad_rect(double u0, double v0, double u1, double v1, double quad[4][2], double QR[3][3])
{
    int ret;
    double du, dv;
    double RQ[3][3];            /* rect->quad transform */

    du = u1-u0;
    dv = v1-v0;
    if (du==0. || dv==0.) {
        //ERROR("pmap_quad_rect: null rectangle\n");
        return PMAP_BAD;
    }

    /* first find mapping from unit uv square to xy quadrilateral */
    ret = pmap_square_quad(quad, RQ);

    qDebug() << RQ[0][0] << RQ[0][1] << RQ[0][2];
    qDebug() << RQ[1][0] << RQ[1][1] << RQ[1][2];
    qDebug() << RQ[2][0] << RQ[2][1] << RQ[2][2];

    if (ret==PMAP_BAD) return PMAP_BAD;

    /* concatenate transform from uv rectangle (u0,v0,u1,v1) to unit square */
    RQ[0][0] /= du;
    RQ[1][0] /= dv;
    RQ[2][0] -= RQ[0][0]*u0 + RQ[1][0]*v0;
    RQ[0][1] /= du;
    RQ[1][1] /= dv;
    RQ[2][1] -= RQ[0][1]*u0 + RQ[1][1]*v0;
    RQ[0][2] /= du;
    RQ[1][2] /= dv;
    RQ[2][2] -= RQ[0][2]*u0 + RQ[1][2]*v0;

    qDebug() << RQ[0][0] << RQ[0][1] << RQ[0][2];
    qDebug() << RQ[1][0] << RQ[1][1] << RQ[1][2];
    qDebug() << RQ[2][0] << RQ[2][1] << RQ[2][2];

    /* now RQ is transform from uv rectangle to xy quadrilateral */
    /* QR = inverse transform, which maps xy to uv */
    if (mx3d_adjoint(RQ, QR)==0.)
    {
        qDebug() << "pmap_quad_rect: warning: determinant=0";
        //ERROR("pmap_quad_rect: warning: determinant=0\n");
    }
    return ret;
}

/*
 * pmap_square_quad: find mapping between unit square and quadrilateral.
 * The correspondence is:
 *
 *      (0,0) --> quad[0]
 *      (1,0) --> quad[1]
 *      (1,1) --> quad[2]
 *      (0,1) --> quad[3]
 */

int pmap_square_quad(double quad[4][2], double SQ[3][3])
{
    double px, py;

    px = X(0)-X(1)+X(2)-X(3);
    py = Y(0)-Y(1)+Y(2)-Y(3);

    qDebug() << "px: " << px << "py: " << py;

    if (ZERO(px) && ZERO(py)){         /* affine */
        qDebug() << "affine";
        SQ[0][0] = X(1)-X(0);
        SQ[1][0] = X(2)-X(1);
        SQ[2][0] = X(0);
        SQ[0][1] = Y(1)-Y(0);
        SQ[1][1] = Y(2)-Y(1);
        SQ[2][1] = Y(0);
        SQ[0][2] = 0.;
        SQ[1][2] = 0.;
        SQ[2][2] = 1.;
        return PMAP_AFFINE;
    }
    else {                              /* projective */
        qDebug() << "projective";
        double dx1, dx2, dy1, dy2, del;

        dx1 = X(1)-X(2);
        dx2 = X(3)-X(2);
        dy1 = Y(1)-Y(2);
        dy2 = Y(3)-Y(2);
        del = DET2(dx1,dx2, dy1,dy2);
        if (del==0.) {
            //ERROR("pmap_square_quad: bad mapping\n");
            return PMAP_BAD;
        }
        SQ[0][2] = DET2(px,dx2, py,dy2)/del;
        SQ[1][2] = DET2(dx1,px, dy1,py)/del;
        SQ[2][2] = 1.;
        SQ[0][0] = X(1)-X(0)+SQ[0][2]*X(1);
        SQ[1][0] = X(3)-X(0)+SQ[1][2]*X(3);
        SQ[2][0] = X(0);
        SQ[0][1] = Y(1)-Y(0)+SQ[0][2]*Y(1);
        SQ[1][1] = Y(3)-Y(0)+SQ[1][2]*Y(3);
        SQ[2][1] = Y(0);
        return PMAP_PROJECTIVE;
    }
}

int pmap_poly(poly *p, double ST[3][3])
{
    double scr[4][2];                   /* vertices of screen quadrilateral */

    //if (p->n!=4)
        //ERROR("only do quadrilaterals at the moment\n");

    /* if edges 0-1 and 2-3 are horz, 1-2 and 3-0 are vert */
    if (V(0)==V(1) && V(2)==V(3) && U(1)==U(2) && U(3)==U(0)) {
        qDebug() << " if edges 0-1 and 2-3 are horz, 1-2 and 3-0 are vert";
        scr[0][0] = X(0); scr[0][1] = Y(0);
        scr[1][0] = X(1); scr[1][1] = Y(1);
        scr[2][0] = X(2); scr[2][1] = Y(2);
        scr[3][0] = X(3); scr[3][1] = Y(3);
        return pmap_quad_rect(U(0), V(0), U(2), V(2), scr, ST);
    }

    /* if edges 0-1 and 2-3 are vert, 1-2 and 3-0 are horz */
    else if (U(0)==U(1) && U(2)==U(3) && V(1)==V(2) && V(3)==V(0)) {
        qDebug() << "if edges 0-1 and 2-3 are vert, 1-2 and 3-0 are horz";
        scr[0][0] = X(1); scr[0][1] = Y(1);
        scr[1][0] = X(2); scr[1][1] = Y(2);
        scr[2][0] = X(3); scr[2][1] = Y(3);
        scr[3][0] = X(0); scr[3][1] = Y(0);
        return pmap_quad_rect(U(1), V(1), U(3), V(3), scr, ST);
    }

    /* if texture is not an orthogonally-oriented rectangle */
    else
    {
        qDebug() << "if texture is not an orthogonally-oriented rectangle";
        return pmap_quad_quad(p, ST);
    }
}

/*
 * pmap_quad_quad: find the projective mapping ST from screen to texture space
 * given the screen and texture coordinates at the vertices of a quadrilateral.
 *
 * Method: concatenate screen_quad->mid_square and
 * mid_square->texture_quad transform matrices.
 * The alternate method, solving an 8x8 system of equations, takes longer.
 */

int pmap_quad_quad(poly *p, double ST[3][3])
{
    int type1, type2;
    double quad[4][2], MS[3][3];
    double SM[3][3], MT[3][3];          /* screen->mid and mid->texture */

    quad[0][0] = X(0); quad[0][1] = Y(0);
    quad[1][0] = X(1); quad[1][1] = Y(1);
    quad[2][0] = X(2); quad[2][1] = Y(2);
    quad[3][0] = X(3); quad[3][1] = Y(3);
    type1 = pmap_square_quad(quad, MS);
    if (mx3d_adjoint(MS, SM) == 0.)
    {
        //ERROR("pmap_quad_quad: warning: determinant=0\n");
    }

    quad[0][0] = U(0); quad[0][1] = V(0);
    quad[1][0] = U(1); quad[1][1] = V(1);
    quad[2][0] = U(2); quad[2][1] = V(2);
    quad[3][0] = U(3); quad[3][1] = V(3);
    type2 = pmap_square_quad(quad, MT);

    if (type1==PMAP_BAD || type2==PMAP_BAD) return PMAP_BAD;

    mx3d_mul(SM, MT, ST);

    return type1|type2;                 /* PMAP_PROJECTIVE prevails */
}

double mx3d_adjoint(double a[3][3], double b[3][3])
{
    b[0][0] = DET2(a[1][1], a[1][2], a[2][1], a[2][2]);
    b[1][0] = DET2(a[1][2], a[1][0], a[2][2], a[2][0]);
    b[2][0] = DET2(a[1][0], a[1][1], a[2][0], a[2][1]);
    b[0][1] = DET2(a[2][1], a[2][2], a[0][1], a[0][2]);
    b[1][1] = DET2(a[2][2], a[2][0], a[0][2], a[0][0]);
    b[2][1] = DET2(a[2][0], a[2][1], a[0][0], a[0][1]);
    b[0][2] = DET2(a[0][1], a[0][2], a[1][1], a[1][2]);
    b[1][2] = DET2(a[0][2], a[0][0], a[1][2], a[1][0]);
    b[2][2] = DET2(a[0][0], a[0][1], a[1][0], a[1][1]);

    qDebug() << b[0][0] << b[0][1] << b[0][2];
    qDebug() << b[1][0] << b[1][1] << b[1][2];
    qDebug() << b[2][0] << b[2][1] << b[2][2];

    return a[0][0]*b[0][0] + a[0][1]*b[0][1] + a[0][2]*b[0][2];
}

/* mx3_mul: matrix multiply: c = a*b */

void mx3d_mul(double a[3][3], double b[3][3], double c[3][3])
{
    c[0][0] = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0];
    c[0][1] = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1];
    c[0][2] = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2];
    c[1][0] = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0];
    c[1][1] = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1];
    c[1][2] = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2];
    c[2][0] = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0];
    c[2][1] = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1];
    c[2][2] = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
}

/* mx3d_transform: transform point p by matrix a: q = p*a */

void mx3d_transform(double p[3], double a[3][3], double q[3])
{
    q[0] = p[0]*a[0][0] + p[1]*a[1][0] + p[2]*a[2][0];
    q[1] = p[0]*a[0][1] + p[1]*a[1][1] + p[2]*a[2][1];
    q[2] = p[0]*a[0][2] + p[1]*a[1][2] + p[2]*a[2][2];
}