src/px/pmapgen.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 * pmapgen.c: general routines for 2-D perspective mappings.
 * These routines are independent of the poly structure,
 * so we do not think in terms of texture and screen space.
 *
 * Paul Heckbert        5 Nov 85, 12 Dec 85
 */

static char rcsid[] = "$Header: pmapgen.c,v 2.0 88/10/12 21:58:33 ph Locked $";
#include <stdio.h>
#include "pmap.h"
#include "mx3.h"

#define TOLERANCE 1e-13
#define ZERO(x) ((x)<TOLERANCE && (x)>-TOLERANCE)

#define X(i) qdrl[i][0]         /* quadrilateral x and y */
#define Y(i) qdrl[i][1]

/*
 * pmap_quad_rect: find mapping between quadrilateral and rectangle.
 * The correspondence is:
 *
 *      qdrl[0] --> (u0,v0)
 *      qdrl[1] --> (u1,v0)
 *      qdrl[2] --> (u1,v1)
 *      qdrl[3] --> (u0,v1)
 *
 * This method of computing the adjoint numerically is cheaper than
 * computing it symbolically.
 */
        
pmap_quad_rect(u0, v0, u1, v1, qdrl, QR)
double u0, v0, u1, v1;          /* bounds of rectangle */
double qdrl[4][2];              /* vertices of quadrilateral */
double QR[3][3];                /* qdrl->rect transform (returned) */
{
    int ret;
    double du, dv;
    double RQ[3][3];            /* rect->qdrl transform */

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

    /* first find mapping from unit uv square to xy quadrilateral */
    ret = pmap_square_quad(qdrl, RQ);
    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;

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

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

pmap_square_quad(qdrl, SQ)
register double qdrl[4][2];     /* vertices of quadrilateral */
register double SQ[3][3];       /* square->qdrl transform */
{
    double px, py;

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

    if (ZERO(px) && ZERO(py)) {         /* 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_LINEAR;
    }
    else {                              /* perspective */
        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.) {
            fprintf(stderr, "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_PERSP;
    }
}
Revision:	2.3
Committed:	Sat Feb 22 02:07:27 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R5
Changes since 2.2:	+1 -4 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.3	static const char RCSid[] = "$Id$";
3	greg	2.1	#endif
4			/*
5			* pmapgen.c: general routines for 2-D perspective mappings.
6			* These routines are independent of the poly structure,
7			* so we do not think in terms of texture and screen space.
8			*
9			* Paul Heckbert 5 Nov 85, 12 Dec 85
10			*/
11
12			static char rcsid[] = "$Header: pmapgen.c,v 2.0 88/10/12 21:58:33 ph Locked $";
13			#include <stdio.h>
14			#include "pmap.h"
15			#include "mx3.h"
16
17			#define TOLERANCE 1e-13
18			#define ZERO(x) ((x)<TOLERANCE && (x)>-TOLERANCE)
19
20	greg	2.2	#define X(i) qdrl[i][0] /* quadrilateral x and y */
21			#define Y(i) qdrl[i][1]
22	greg	2.1
23			/*
24			* pmap_quad_rect: find mapping between quadrilateral and rectangle.
25			* The correspondence is:
26			*
27	greg	2.2	* qdrl[0] --> (u0,v0)
28			* qdrl[1] --> (u1,v0)
29			* qdrl[2] --> (u1,v1)
30			* qdrl[3] --> (u0,v1)
31	greg	2.1	*
32			* This method of computing the adjoint numerically is cheaper than
33			* computing it symbolically.
34			*/
35
36	greg	2.2	pmap_quad_rect(u0, v0, u1, v1, qdrl, QR)
37	greg	2.1	double u0, v0, u1, v1; /* bounds of rectangle */
38	greg	2.2	double qdrl[4][2]; /* vertices of quadrilateral */
39			double QR[3][3]; /* qdrl->rect transform (returned) */
40	greg	2.1	{
41			int ret;
42			double du, dv;
43	greg	2.2	double RQ[3][3]; /* rect->qdrl transform */
44	greg	2.1
45			du = u1-u0;
46			dv = v1-v0;
47			if (du==0. \|\| dv==0.) {
48			fprintf(stderr, "pmap_quad_rect: null rectangle\n");
49			return PMAP_BAD;
50			}
51
52			/* first find mapping from unit uv square to xy quadrilateral */
53	greg	2.2	ret = pmap_square_quad(qdrl, RQ);
54	greg	2.1	if (ret==PMAP_BAD) return PMAP_BAD;
55
56			/* concatenate transform from uv rectangle (u0,v0,u1,v1) to unit square */
57			RQ[0][0] /= du;
58			RQ[1][0] /= dv;
59			RQ[2][0] -= RQ[0][0]u0 + RQ[1][0]v0;
60			RQ[0][1] /= du;
61			RQ[1][1] /= dv;
62			RQ[2][1] -= RQ[0][1]u0 + RQ[1][1]v0;
63			RQ[0][2] /= du;
64			RQ[1][2] /= dv;
65			RQ[2][2] -= RQ[0][2]u0 + RQ[1][2]v0;
66
67			/* now RQ is transform from uv rectangle to xy quadrilateral */
68			/* QR = inverse transform, which maps xy to uv */
69			if (mx3d_adjoint(RQ, QR)==0.)
70			fprintf(stderr, "pmap_quad_rect: warning: determinant=0\n");
71			return ret;
72			}
73
74			/*
75			* pmap_square_quad: find mapping between unit square and quadrilateral.
76			* The correspondence is:
77			*
78	greg	2.2	* (0,0) --> qdrl[0]
79			* (1,0) --> qdrl[1]
80			* (1,1) --> qdrl[2]
81			* (0,1) --> qdrl[3]
82	greg	2.1	*/
83
84	greg	2.2	pmap_square_quad(qdrl, SQ)
85			register double qdrl[4][2]; /* vertices of quadrilateral */
86			register double SQ[3][3]; /* square->qdrl transform */
87	greg	2.1	{
88			double px, py;
89
90			px = X(0)-X(1)+X(2)-X(3);
91			py = Y(0)-Y(1)+Y(2)-Y(3);
92
93			if (ZERO(px) && ZERO(py)) { /* affine */
94			SQ[0][0] = X(1)-X(0);
95			SQ[1][0] = X(2)-X(1);
96			SQ[2][0] = X(0);
97			SQ[0][1] = Y(1)-Y(0);
98			SQ[1][1] = Y(2)-Y(1);
99			SQ[2][1] = Y(0);
100			SQ[0][2] = 0.;
101			SQ[1][2] = 0.;
102			SQ[2][2] = 1.;
103			return PMAP_LINEAR;
104			}
105			else { /* perspective */
106			double dx1, dx2, dy1, dy2, del;
107
108			dx1 = X(1)-X(2);
109			dx2 = X(3)-X(2);
110			dy1 = Y(1)-Y(2);
111			dy2 = Y(3)-Y(2);
112			del = DET2(dx1,dx2, dy1,dy2);
113			if (del==0.) {
114			fprintf(stderr, "pmap_square_quad: bad mapping\n");
115			return PMAP_BAD;
116			}
117			SQ[0][2] = DET2(px,dx2, py,dy2)/del;
118			SQ[1][2] = DET2(dx1,px, dy1,py)/del;
119			SQ[2][2] = 1.;
120			SQ[0][0] = X(1)-X(0)+SQ[0][2]*X(1);
121			SQ[1][0] = X(3)-X(0)+SQ[1][2]*X(3);
122			SQ[2][0] = X(0);
123			SQ[0][1] = Y(1)-Y(0)+SQ[0][2]*Y(1);
124			SQ[1][1] = Y(3)-Y(0)+SQ[1][2]*Y(3);
125			SQ[2][1] = Y(0);
126			return PMAP_PERSP;
127			}
128			}