Back to index

radiance  4R0+20100331
pmapgen.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: pmapgen.c,v 2.4 2004/03/28 20:33:14 schorsch Exp $";
00003 #endif
00004 /*
00005  * pmapgen.c: general routines for 2-D perspective mappings.
00006  * These routines are independent of the poly structure,
00007  * so we do not think in terms of texture and screen space.
00008  *
00009  * Paul Heckbert     5 Nov 85, 12 Dec 85
00010  */
00011 
00012 #include <stdio.h>
00013 #include "pmap.h"
00014 #include "mx3.h"
00015 
00016 #define TOLERANCE 1e-13
00017 #define ZERO(x) ((x)<TOLERANCE && (x)>-TOLERANCE)
00018 
00019 #define X(i) qdrl[i][0]            /* quadrilateral x and y */
00020 #define Y(i) qdrl[i][1]
00021 
00022 /*
00023  * pmap_quad_rect: find mapping between quadrilateral and rectangle.
00024  * The correspondence is:
00025  *
00026  *     qdrl[0] --> (u0,v0)
00027  *     qdrl[1] --> (u1,v0)
00028  *     qdrl[2] --> (u1,v1)
00029  *     qdrl[3] --> (u0,v1)
00030  *
00031  * This method of computing the adjoint numerically is cheaper than
00032  * computing it symbolically.
00033  */
00034        
00035 extern int
00036 pmap_quad_rect(
00037        double u0,           /* bounds of rectangle */
00038        double v0,
00039        double u1,
00040        double v1,
00041        double qdrl[4][2],          /* vertices of quadrilateral */
00042        double QR[3][3]             /* qdrl->rect transform (returned) */
00043 )
00044 {
00045     int ret;
00046     double du, dv;
00047     double RQ[3][3];        /* rect->qdrl transform */
00048 
00049     du = u1-u0;
00050     dv = v1-v0;
00051     if (du==0. || dv==0.) {
00052        fprintf(stderr, "pmap_quad_rect: null rectangle\n");
00053        return PMAP_BAD;
00054     }
00055 
00056     /* first find mapping from unit uv square to xy quadrilateral */
00057     ret = pmap_square_quad(qdrl, RQ);
00058     if (ret==PMAP_BAD) return PMAP_BAD;
00059 
00060     /* concatenate transform from uv rectangle (u0,v0,u1,v1) to unit square */
00061     RQ[0][0] /= du;
00062     RQ[1][0] /= dv;
00063     RQ[2][0] -= RQ[0][0]*u0 + RQ[1][0]*v0;
00064     RQ[0][1] /= du;
00065     RQ[1][1] /= dv;
00066     RQ[2][1] -= RQ[0][1]*u0 + RQ[1][1]*v0;
00067     RQ[0][2] /= du;
00068     RQ[1][2] /= dv;
00069     RQ[2][2] -= RQ[0][2]*u0 + RQ[1][2]*v0;
00070 
00071     /* now RQ is transform from uv rectangle to xy quadrilateral */
00072     /* QR = inverse transform, which maps xy to uv */
00073     if (mx3d_adjoint(RQ, QR)==0.)
00074        fprintf(stderr, "pmap_quad_rect: warning: determinant=0\n");
00075     return ret;
00076 }
00077 
00078 /*
00079  * pmap_square_quad: find mapping between unit square and quadrilateral.
00080  * The correspondence is:
00081  *
00082  *     (0,0) --> qdrl[0]
00083  *     (1,0) --> qdrl[1]
00084  *     (1,1) --> qdrl[2]
00085  *     (0,1) --> qdrl[3]
00086  */
00087 
00088 extern int
00089 pmap_square_quad(
00090        register double qdrl[4][2], /* vertices of quadrilateral */
00091        register double SQ[3][3]    /* square->qdrl transform */
00092 )
00093 {
00094     double px, py;
00095 
00096     px = X(0)-X(1)+X(2)-X(3);
00097     py = Y(0)-Y(1)+Y(2)-Y(3);
00098 
00099     if (ZERO(px) && ZERO(py)) {           /* affine */
00100        SQ[0][0] = X(1)-X(0);
00101        SQ[1][0] = X(2)-X(1);
00102        SQ[2][0] = X(0);
00103        SQ[0][1] = Y(1)-Y(0);
00104        SQ[1][1] = Y(2)-Y(1);
00105        SQ[2][1] = Y(0);
00106        SQ[0][2] = 0.;
00107        SQ[1][2] = 0.;
00108        SQ[2][2] = 1.;
00109        return PMAP_LINEAR;
00110     }
00111     else {                         /* perspective */
00112        double dx1, dx2, dy1, dy2, del;
00113 
00114        dx1 = X(1)-X(2);
00115        dx2 = X(3)-X(2);
00116        dy1 = Y(1)-Y(2);
00117        dy2 = Y(3)-Y(2);
00118        del = DET2(dx1,dx2, dy1,dy2);
00119        if (del==0.) {
00120            fprintf(stderr, "pmap_square_quad: bad mapping\n");
00121            return PMAP_BAD;
00122        }
00123        SQ[0][2] = DET2(px,dx2, py,dy2)/del;
00124        SQ[1][2] = DET2(dx1,px, dy1,py)/del;
00125        SQ[2][2] = 1.;
00126        SQ[0][0] = X(1)-X(0)+SQ[0][2]*X(1);
00127        SQ[1][0] = X(3)-X(0)+SQ[1][2]*X(3);
00128        SQ[2][0] = X(0);
00129        SQ[0][1] = Y(1)-Y(0)+SQ[0][2]*Y(1);
00130        SQ[1][1] = Y(3)-Y(0)+SQ[1][2]*Y(3);
00131        SQ[2][1] = Y(0);
00132        return PMAP_PERSP;
00133     }
00134 }