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
#	Content
1	#ifndef lint
2	static const char RCSid[] = "$Id$";
3	#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	#define X(i) qdrl[i][0] /* quadrilateral x and y */
21	#define Y(i) qdrl[i][1]
22
23	/*
24	* pmap_quad_rect: find mapping between quadrilateral and rectangle.
25	* The correspondence is:
26	*
27	* qdrl[0] --> (u0,v0)
28	* qdrl[1] --> (u1,v0)
29	* qdrl[2] --> (u1,v1)
30	* qdrl[3] --> (u0,v1)
31	*
32	* This method of computing the adjoint numerically is cheaper than
33	* computing it symbolically.
34	*/
35
36	pmap_quad_rect(u0, v0, u1, v1, qdrl, QR)
37	double u0, v0, u1, v1; /* bounds of rectangle */
38	double qdrl[4][2]; /* vertices of quadrilateral */
39	double QR[3][3]; /* qdrl->rect transform (returned) */
40	{
41	int ret;
42	double du, dv;
43	double RQ[3][3]; /* rect->qdrl transform */
44
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	ret = pmap_square_quad(qdrl, RQ);
54	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	* (0,0) --> qdrl[0]
79	* (1,0) --> qdrl[1]
80	* (1,1) --> qdrl[2]
81	* (0,1) --> qdrl[3]
82	*/
83
84	pmap_square_quad(qdrl, SQ)
85	register double qdrl[4][2]; /* vertices of quadrilateral */
86	register double SQ[3][3]; /* square->qdrl transform */
87	{
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	}