Home | History | Annotate | Line # | Download | only in internals
      1 /*
      2 ** License Applicability. Except to the extent portions of this file are
      3 ** made subject to an alternative license as permitted in the SGI Free
      4 ** Software License B, Version 1.1 (the "License"), the contents of this
      5 ** file are subject only to the provisions of the License. You may not use
      6 ** this file except in compliance with the License. You may obtain a copy
      7 ** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
      8 ** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
      9 **
     10 ** http://oss.sgi.com/projects/FreeB
     11 **
     12 ** Note that, as provided in the License, the Software is distributed on an
     13 ** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
     14 ** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
     15 ** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
     16 ** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
     17 **
     18 ** Original Code. The Original Code is: OpenGL Sample Implementation,
     19 ** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
     20 ** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
     21 ** Copyright in any portions created by third parties is as indicated
     22 ** elsewhere herein. All Rights Reserved.
     23 **
     24 ** Additional Notice Provisions: The application programming interfaces
     25 ** established by SGI in conjunction with the Original Code are The
     26 ** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
     27 ** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
     28 ** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
     29 ** Window System(R) (Version 1.3), released October 19, 1998. This software
     30 ** was created using the OpenGL(R) version 1.2.1 Sample Implementation
     31 ** published by SGI, but has not been independently verified as being
     32 ** compliant with the OpenGL(R) version 1.2.1 Specification.
     33 */
     34 
     35 /*
     36  * mapdesc.c++
     37  *
     38  */
     39 
     40 #include <stdio.h>
     41 #include "glimports.h"
     42 #include "mystdio.h"
     43 #include "myassert.h"
     44 #include "mystring.h"
     45 #include "mymath.h"
     46 #include "backend.h"
     47 #include "nurbsconsts.h"
     48 #include "mapdesc.h"
     49 
     50 Mapdesc::Mapdesc( long _type, int _israt, int _ncoords, Backend& b )
     51     : backend( b )
     52 {
     53     type 		= _type;
     54     isrational 		= _israt;
     55     ncoords 		= _ncoords;
     56     hcoords		= _ncoords + (_israt ? 0 : 1 );
     57     inhcoords		= _ncoords - (_israt ? 1 : 0 );
     58     mask 		= ((1<<(inhcoords*2))-1);
     59     next		= 0;
     60 
     61     assert( hcoords <= MAXCOORDS );
     62     assert( inhcoords >= 1 );
     63 
     64     pixel_tolerance 	= 1.0;
     65     error_tolerance	= 1.0;
     66     bbox_subdividing	= N_NOBBOXSUBDIVISION;
     67     culling_method 	= N_NOCULLING;
     68     sampling_method 	= N_NOSAMPLING;
     69     clampfactor 	= N_NOCLAMPING;
     70     minsavings 		= N_NOSAVINGSSUBDIVISION;
     71     s_steps  		= 0.0;
     72     t_steps 		= 0.0;
     73     maxrate 		= ( s_steps < 0.0 ) ? 0.0 : s_steps;
     74     maxsrate 		= ( s_steps < 0.0 ) ? 0.0 : s_steps;
     75     maxtrate 		= ( t_steps < 0.0 ) ? 0.0 : t_steps;
     76     identify( bmat );
     77     identify( cmat );
     78     identify( smat );
     79     for( int i = 0; i != inhcoords; i++ )
     80 	bboxsize[i] = 1.0;
     81 }
     82 
     83 void
     84 Mapdesc::setBboxsize( INREAL *mat )
     85 {
     86     for( int i = 0; i != inhcoords; i++ )
     87 	bboxsize[i] = (REAL) mat[i];
     88 }
     89 
     90 void
     91 Mapdesc::identify( REAL dest[MAXCOORDS][MAXCOORDS] )
     92 {
     93     memset( dest, 0, sizeof( REAL ) * MAXCOORDS * MAXCOORDS );
     94     for( int i=0; i != hcoords; i++ )
     95 	dest[i][i] = 1.0;
     96 }
     97 
     98 void
     99 Mapdesc::surfbbox( REAL bb[2][MAXCOORDS] )
    100 {
    101     backend.surfbbox( type, bb[0], bb[1] );
    102 }
    103 
    104 void
    105 Mapdesc::copy( REAL dest[MAXCOORDS][MAXCOORDS], long n, INREAL *src,
    106 	long rstride, long cstride )
    107 {
    108     assert( n >= 0 );
    109     for( int i=0; i != n; i++ )
    110         for( int j=0; j != n; j++ )
    111 	    dest[i][j] = src[i*rstride + j*cstride];
    112 }
    113 
    114 /*--------------------------------------------------------------------------
    115  * copyPt - copy a homogeneous point
    116  *--------------------------------------------------------------------------
    117  */
    118 void
    119 Mapdesc::copyPt( REAL *d, REAL *s )
    120 {
    121     assert( hcoords > 0 );
    122     switch( hcoords  ) {
    123 	case 4:
    124 	    d[3] = s[3];
    125 	    d[2] = s[2];
    126 	    d[1] = s[1];
    127 	    d[0] = s[0];
    128 	    break;
    129 	case 3:
    130 	    d[2] = s[2];
    131 	    d[1] = s[1];
    132 	    d[0] = s[0];
    133 	    break;
    134 	case 2:
    135 	    d[1] = s[1];
    136 	    d[0] = s[0];
    137 	    break;
    138 	case 1:
    139 	    d[0] = s[0];
    140 	    break;
    141 	case 5:
    142 	    d[4] = s[4];
    143 	    d[3] = s[3];
    144 	    d[2] = s[2];
    145 	    d[1] = s[1];
    146 	    d[0] = s[0];
    147 	    break;
    148 	default:
    149 	    memcpy( d, s, hcoords * sizeof( REAL ) );
    150 	    break;
    151     }
    152 }
    153 
    154 /*--------------------------------------------------------------------------
    155  * sumPt - compute affine combination of two homogeneous points
    156  *--------------------------------------------------------------------------
    157  */
    158 void
    159 Mapdesc::sumPt( REAL *dst, REAL *src1, REAL *src2, REAL alpha, REAL beta )
    160 {
    161     assert( hcoords > 0 );
    162     switch( hcoords  ) {
    163 	case 4:
    164 	    dst[3] = src1[3] * alpha + src2[3] * beta;
    165 	    dst[2] = src1[2] * alpha + src2[2] * beta;
    166 	    dst[1] = src1[1] * alpha + src2[1] * beta;
    167 	    dst[0] = src1[0] * alpha + src2[0] * beta;
    168 	    break;
    169 	case 3:
    170 	    dst[2] = src1[2] * alpha + src2[2] * beta;
    171 	    dst[1] = src1[1] * alpha + src2[1] * beta;
    172 	    dst[0] = src1[0] * alpha + src2[0] * beta;
    173 	    break;
    174 	case 2:
    175 	    dst[1] = src1[1] * alpha + src2[1] * beta;
    176 	    dst[0] = src1[0] * alpha + src2[0] * beta;
    177 	    break;
    178 	case 1:
    179 	    dst[0] = src1[0] * alpha + src2[0] * beta;
    180 	    break;
    181 	case 5:
    182 	    dst[4] = src1[4] * alpha + src2[4] * beta;
    183 	    dst[3] = src1[3] * alpha + src2[3] * beta;
    184 	    dst[2] = src1[2] * alpha + src2[2] * beta;
    185 	    dst[1] = src1[1] * alpha + src2[1] * beta;
    186 	    dst[0] = src1[0] * alpha + src2[0] * beta;
    187 	    break;
    188 	default: {
    189 		for( int i = 0; i != hcoords; i++ )
    190 		    dst[i] = src1[i] * alpha + src2[i] * beta;
    191             }
    192 	    break;
    193     }
    194 }
    195 
    196 /*--------------------------------------------------------------------------
    197  * clipbits - compute bit-vector indicating point/window position
    198  *		       of a (transformed) homogeneous point
    199  *--------------------------------------------------------------------------
    200  */
    201 unsigned int
    202 Mapdesc::clipbits( REAL *p )
    203 {
    204     assert( inhcoords >= 0 );
    205     assert( inhcoords <= 3 );
    206 
    207     int nc = inhcoords;
    208     REAL pw = p[nc];
    209     REAL nw = -pw;
    210     unsigned int bits = 0;
    211 
    212     if( pw == 0.0 ) return mask;
    213 
    214     if( pw > 0.0 ) {
    215 	switch( nc ) {
    216 	case 3:
    217 	    if( p[2] <= pw ) bits |= (1<<5);
    218 	    if( p[2] >= nw ) bits |= (1<<4);
    219 	    if( p[1] <= pw ) bits |= (1<<3);
    220 	    if( p[1] >= nw ) bits |= (1<<2);
    221 	    if( p[0] <= pw ) bits |= (1<<1);
    222 	    if( p[0] >= nw ) bits |= (1<<0);
    223             return bits;
    224 	case 2:
    225 	    if( p[1] <= pw ) bits |= (1<<3);
    226 	    if( p[1] >= nw ) bits |= (1<<2);
    227 	    if( p[0] <= pw ) bits |= (1<<1);
    228 	    if( p[0] >= nw ) bits |= (1<<0);
    229             return bits;
    230 	case 1:
    231 	    if( p[0] <= pw ) bits |= (1<<1);
    232 	    if( p[0] >= nw ) bits |= (1<<0);
    233             return bits;
    234 	default: {
    235 		int bit = 1;
    236 		for( int i=0; i<nc; i++ ) {
    237 		    if( p[i] >= nw ) bits |= bit;
    238 		    bit <<= 1;
    239 		    if( p[i] <= pw ) bits |= bit;
    240 		    bit <<= 1;
    241 		}
    242 		abort();
    243 		break;
    244 	    }
    245 	}
    246     } else {
    247 	switch( nc ) {
    248 	case 3:
    249 	    if( p[2] <= nw ) bits |= (1<<5);
    250 	    if( p[2] >= pw ) bits |= (1<<4);
    251 	    if( p[1] <= nw ) bits |= (1<<3);
    252 	    if( p[1] >= pw ) bits |= (1<<2);
    253 	    if( p[0] <= nw ) bits |= (1<<1);
    254 	    if( p[0] >= pw ) bits |= (1<<0);
    255             return bits;
    256 	case 2:
    257 	    if( p[1] <= nw ) bits |= (1<<3);
    258 	    if( p[1] >= pw ) bits |= (1<<2);
    259 	    if( p[0] <= nw ) bits |= (1<<1);
    260 	    if( p[0] >= pw ) bits |= (1<<0);
    261             return bits;
    262 	case 1:
    263 	    if( p[0] <= nw ) bits |= (1<<1);
    264 	    if( p[0] >= pw ) bits |= (1<<0);
    265             return bits;
    266 	default: {
    267 		int bit = 1;
    268 		for( int i=0; i<nc; i++ ) {
    269 		    if( p[i] >= pw ) bits |= bit;
    270 		    bit <<= 1;
    271 		    if( p[i] <= nw ) bits |= bit;
    272 		    bit <<= 1;
    273 		}
    274 		abort();
    275 		break;
    276 	    }
    277 	}
    278     }
    279     return bits;
    280 }
    281 
    282 /*--------------------------------------------------------------------------
    283  * xformRational - transform a homogeneous point
    284  *--------------------------------------------------------------------------
    285  */
    286 void
    287 Mapdesc::xformRational( Maxmatrix mat, REAL *d, REAL *s )
    288 {
    289     assert( hcoords >= 0 );
    290 
    291     if( hcoords == 3 ) {
    292 	REAL x = s[0];
    293 	REAL y = s[1];
    294 	REAL z = s[2];
    295 	d[0] = x*mat[0][0]+y*mat[1][0]+z*mat[2][0];
    296 	d[1] = x*mat[0][1]+y*mat[1][1]+z*mat[2][1];
    297 	d[2] = x*mat[0][2]+y*mat[1][2]+z*mat[2][2];
    298     } else if( hcoords == 4 ) {
    299 	REAL x = s[0];
    300 	REAL y = s[1];
    301 	REAL z = s[2];
    302 	REAL w = s[3];
    303 	d[0] = x*mat[0][0]+y*mat[1][0]+z*mat[2][0]+w*mat[3][0];
    304 	d[1] = x*mat[0][1]+y*mat[1][1]+z*mat[2][1]+w*mat[3][1];
    305 	d[2] = x*mat[0][2]+y*mat[1][2]+z*mat[2][2]+w*mat[3][2];
    306 	d[3] = x*mat[0][3]+y*mat[1][3]+z*mat[2][3]+w*mat[3][3];
    307     } else {
    308 	for( int i=0; i != hcoords; i++ ) {
    309 	    d[i] = 0;
    310 	    for( int j = 0; j != hcoords; j++ )
    311 		d[i] += s[j] * mat[j][i];
    312 	}
    313     }
    314 }
    315 
    316 /*--------------------------------------------------------------------------
    317  * xformNonrational - transform a inhomogeneous point to a homogeneous point
    318  *--------------------------------------------------------------------------
    319  */
    320 void
    321 Mapdesc::xformNonrational( Maxmatrix mat, REAL *d, REAL *s )
    322 {
    323     if( inhcoords == 2 ) {
    324 	REAL x = s[0];
    325 	REAL y = s[1];
    326 	d[0] = x*mat[0][0]+y*mat[1][0]+mat[2][0];
    327 	d[1] = x*mat[0][1]+y*mat[1][1]+mat[2][1];
    328 	d[2] = x*mat[0][2]+y*mat[1][2]+mat[2][2];
    329     } else if( inhcoords == 3 ) {
    330 	REAL x = s[0];
    331 	REAL y = s[1];
    332 	REAL z = s[2];
    333 	d[0] = x*mat[0][0]+y*mat[1][0]+z*mat[2][0]+mat[3][0];
    334 	d[1] = x*mat[0][1]+y*mat[1][1]+z*mat[2][1]+mat[3][1];
    335 	d[2] = x*mat[0][2]+y*mat[1][2]+z*mat[2][2]+mat[3][2];
    336 	d[3] = x*mat[0][3]+y*mat[1][3]+z*mat[2][3]+mat[3][3];
    337     } else {
    338         assert( inhcoords >= 0 );
    339 	for( int i=0; i != hcoords; i++ ) {
    340 	    d[i] = mat[inhcoords][i];
    341 	    for( int j = 0; j < inhcoords; j++ )
    342 		d[i] += s[j] * mat[j][i];
    343 	}
    344     }
    345 }
    346 
    347 /*--------------------------------------------------------------------------
    348  * xformAndCullCheck - transform a set of points that may be EITHER
    349  *	homogeneous or inhomogeneous depending on the map description and
    350  *	check if they are either completely inside, completely outside,
    351  *	or intersecting the viewing frustrum.
    352  *--------------------------------------------------------------------------
    353  */
    354 int
    355 Mapdesc::xformAndCullCheck(
    356     REAL *pts, int uorder, int ustride, int vorder, int vstride )
    357 {
    358     assert( uorder > 0 );
    359     assert( vorder > 0 );
    360 
    361     unsigned int inbits = mask;
    362     unsigned int outbits = 0;
    363 
    364     REAL *p = pts;
    365     for( REAL *pend = p + uorder * ustride; p != pend; p += ustride ) {
    366         REAL *q = p;
    367         for( REAL *qend = q + vorder * vstride; q != qend; q += vstride ) {
    368     	    REAL cpts[MAXCOORDS];
    369 	    xformCulling( cpts, q );
    370 	    unsigned int bits = clipbits( cpts );
    371 	    outbits |= bits;
    372 	    inbits &= bits;
    373 	    if( ( outbits == (unsigned int)mask ) && ( inbits != (unsigned int)mask ) ) return CULL_ACCEPT;
    374 	}
    375     }
    376 
    377     if( outbits != (unsigned int)mask ) {
    378 	return CULL_TRIVIAL_REJECT;
    379     } else if( inbits == (unsigned int)mask ) {
    380 	return CULL_TRIVIAL_ACCEPT;
    381     } else {
    382 	return CULL_ACCEPT;
    383     }
    384 }
    385 
    386 /*--------------------------------------------------------------------------
    387  * cullCheck - check if a set of homogeneous transformed points are
    388  *	either completely inside, completely outside,
    389  *	or intersecting the viewing frustrum.
    390  *--------------------------------------------------------------------------
    391  */
    392 int
    393 Mapdesc::cullCheck( REAL *pts, int uorder, int ustride, int vorder, int vstride )
    394 {
    395     unsigned int inbits = mask;
    396     unsigned int outbits  = 0;
    397 
    398     REAL *p = pts;
    399     for( REAL *pend = p + uorder * ustride; p != pend; p += ustride ) {
    400         REAL *q = p;
    401         for( REAL *qend = q + vorder * vstride; q != qend; q += vstride ) {
    402 	    unsigned int bits = clipbits( q );
    403 	    outbits |= bits;
    404 	    inbits &= bits;
    405 	    if( ( outbits == (unsigned int)mask ) && ( inbits != (unsigned int)mask ) ) return CULL_ACCEPT;
    406 	}
    407     }
    408 
    409     if( outbits != (unsigned int)mask ) {
    410 	return CULL_TRIVIAL_REJECT;
    411     } else if( inbits == (unsigned int)mask ) {
    412 	return CULL_TRIVIAL_ACCEPT;
    413     } else {
    414 	return CULL_ACCEPT;
    415     }
    416 }
    417 
    418 /*--------------------------------------------------------------------------
    419  * cullCheck - check if a set of homogeneous transformed points are
    420  *	either completely inside, completely outside,
    421  *	or intersecting the viewing frustrum.
    422  *--------------------------------------------------------------------------
    423  */
    424 int
    425 Mapdesc::cullCheck( REAL *pts, int order, int stride )
    426 {
    427     unsigned int inbits = mask;
    428     unsigned int outbits  = 0;
    429 
    430     REAL *p = pts;
    431     for( REAL *pend = p + order * stride; p != pend; p += stride ) {
    432 	unsigned int bits = clipbits( p );
    433 	outbits |= bits;
    434 	inbits &= bits;
    435 	if( ( outbits == (unsigned int)mask ) && ( inbits != (unsigned int)mask ) ) return CULL_ACCEPT;
    436     }
    437 
    438     if( outbits != (unsigned int)mask ) {
    439 	return CULL_TRIVIAL_REJECT;
    440     } else if( inbits == (unsigned int)mask ) {
    441 	return CULL_TRIVIAL_ACCEPT;
    442     } else {
    443 	return CULL_ACCEPT;
    444     }
    445 }
    446 
    447 /*--------------------------------------------------------------------------
    448  * xformSampling - transform a set of points that may be EITHER
    449  *	homogeneous or inhomogeneous depending on the map description
    450  *	into sampling space
    451  *--------------------------------------------------------------------------
    452  */
    453 void
    454 Mapdesc::xformSampling( REAL *pts, int order, int stride, REAL *sp, int outstride )
    455 {
    456     xformMat( smat, pts, order, stride, sp, outstride );
    457 }
    458 
    459 void
    460 Mapdesc::xformBounding( REAL *pts, int order, int stride, REAL *sp, int outstride )
    461 {
    462     xformMat( bmat, pts, order, stride, sp, outstride );
    463 }
    464 
    465 /*--------------------------------------------------------------------------
    466  * xformCulling - transform a set of points that may be EITHER
    467  *	homogeneous or inhomogeneous depending on the map description
    468  *	into culling space
    469  *--------------------------------------------------------------------------
    470  */
    471 void
    472 Mapdesc::xformCulling( REAL *pts, int order, int stride, REAL *cp, int outstride )
    473 {
    474     xformMat( cmat, pts, order, stride, cp, outstride );
    475 }
    476 
    477 /*--------------------------------------------------------------------------
    478  * xformCulling - transform a set of points that may be EITHER
    479  *	homogeneous or inhomogeneous depending on the map description
    480  *	into culling space
    481  *--------------------------------------------------------------------------
    482  */
    483 void
    484 Mapdesc::xformCulling( REAL *pts,
    485     int uorder, int ustride,
    486     int vorder, int vstride,
    487     REAL *cp, int outustride, int outvstride )
    488 {
    489     xformMat( cmat, pts, uorder, ustride, vorder, vstride, cp, outustride, outvstride );
    490 }
    491 
    492 /*--------------------------------------------------------------------------
    493  * xformSampling - transform a set of points that may be EITHER
    494  *	homogeneous or inhomogeneous depending on the map description
    495  *	into sampling space
    496  *--------------------------------------------------------------------------
    497  */
    498 void
    499 Mapdesc::xformSampling( REAL *pts,
    500     int uorder, int ustride,
    501     int vorder, int vstride,
    502     REAL *sp, int outustride, int outvstride )
    503 {
    504     xformMat( smat, pts, uorder, ustride, vorder, vstride, sp, outustride, outvstride );
    505 }
    506 
    507 void
    508 Mapdesc::xformBounding( REAL *pts,
    509     int uorder, int ustride,
    510     int vorder, int vstride,
    511     REAL *sp, int outustride, int outvstride )
    512 {
    513     xformMat( bmat, pts, uorder, ustride, vorder, vstride, sp, outustride, outvstride );
    514 }
    515 
    516 void
    517 Mapdesc::xformMat(
    518     Maxmatrix	mat,
    519     REAL *	pts,
    520     int 	order,
    521     int 	stride,
    522     REAL *	cp,
    523     int 	outstride )
    524 {
    525     if( isrational ) {
    526 	REAL *pend = pts + order * stride;
    527 	for( REAL *p = pts ; p != pend; p += stride ) {
    528 	    xformRational( mat, cp, p );
    529 	    cp += outstride;
    530 	}
    531     } else {
    532 	REAL *pend = pts + order * stride;
    533 	for( REAL *p = pts ; p != pend; p += stride ) {
    534 	    xformNonrational( mat, cp, p );
    535 	    cp += outstride;
    536 	}
    537     }
    538 }
    539 
    540 void
    541 Mapdesc::xformMat( Maxmatrix mat, REAL *pts,
    542     int uorder, int ustride,
    543     int vorder, int vstride,
    544     REAL *cp, int outustride, int outvstride )
    545 {
    546     if( isrational ) {
    547 	REAL *pend = pts + uorder * ustride;
    548 	for( REAL *p = pts ; p != pend; p += ustride ) {
    549 	    REAL *cpts2 = cp;
    550 	    REAL *qend = p + vorder * vstride;
    551 	    for( REAL *q = p; q != qend; q += vstride ) {
    552 		xformRational( mat, cpts2, q );
    553 		cpts2 += outvstride;
    554 	    }
    555 	    cp += outustride;
    556 	}
    557     } else {
    558 	REAL *pend = pts + uorder * ustride;
    559 	for( REAL *p = pts ; p != pend; p += ustride ) {
    560 	    REAL *cpts2 = cp;
    561 	    REAL *qend = p + vorder * vstride;
    562 	    for( REAL *q = p; q != qend; q += vstride ) {
    563 		xformNonrational( mat, cpts2, q );
    564 		cpts2 += outvstride;
    565 	    }
    566 	    cp += outustride;
    567 	}
    568     }
    569 }
    570 
    571 /*--------------------------------------------------------------------------
    572  * subdivide - subdivide a curve along an isoparametric line
    573  *--------------------------------------------------------------------------
    574  */
    575 
    576 void
    577 Mapdesc::subdivide( REAL *src, REAL *dst, REAL v, int stride, int order )
    578 {
    579     REAL mv = 1.0 - v;
    580 
    581     for( REAL *send=src+stride*order; src!=send; send-=stride, dst+=stride ) {
    582 	copyPt( dst, src );
    583 	REAL *qpnt = src + stride;
    584 	for( REAL *qp=src; qpnt!=send; qp=qpnt, qpnt+=stride )
    585 	    sumPt( qp, qp, qpnt, mv, v );
    586     }
    587 }
    588 
    589 /*--------------------------------------------------------------------------
    590  * subdivide - subdivide a patch along an isoparametric line
    591  *--------------------------------------------------------------------------
    592  */
    593 
    594 void
    595 Mapdesc::subdivide( REAL *src, REAL *dst, REAL v,
    596     int so, int ss, int to, int ts  )
    597 {
    598     REAL mv = 1.0 - v;
    599 
    600     for( REAL *slast = src+ss*so; src != slast; src += ss, dst += ss ) {
    601 	REAL *sp = src;
    602 	REAL *dp = dst;
    603         for( REAL *send = src+ts*to; sp != send; send -= ts, dp += ts ) {
    604 	    copyPt( dp, sp );
    605 	    REAL *qp = sp;
    606 	    for( REAL *qpnt = sp+ts; qpnt != send; qp = qpnt, qpnt += ts )
    607 	        sumPt( qp, qp, qpnt, mv, v );
    608 	}
    609     }
    610 }
    611 
    612 
    613 #define sign(x)	((x > 0) ? 1 : ((x < 0.0) ? -1 : 0))
    614 
    615 /*--------------------------------------------------------------------------
    616  * project - project a set of homogeneous coordinates into inhomogeneous ones
    617  *--------------------------------------------------------------------------
    618  */
    619 int
    620 Mapdesc::project( REAL *src, int rstride, int cstride,
    621 	          REAL *dest, int trstride, int tcstride,
    622 		  int nrows, int ncols )
    623 {
    624     int s = sign( src[inhcoords] );
    625     REAL *rlast = src + nrows * rstride;
    626     REAL *trptr = dest;
    627     for( REAL *rptr=src; rptr != rlast; rptr+=rstride, trptr+=trstride ) {
    628 	REAL *clast = rptr + ncols * cstride;
    629 	REAL *tcptr = trptr;
    630 	for( REAL *cptr = rptr; cptr != clast; cptr+=cstride, tcptr+=tcstride ) {
    631 	    REAL *coordlast = cptr + inhcoords;
    632 	    if( sign( *coordlast ) != s ) return 0;
    633 	    REAL *tcoord = tcptr;
    634 	    for( REAL *coord = cptr; coord != coordlast; coord++, tcoord++ ) {
    635 		*tcoord = *coord / *coordlast;
    636 	    }
    637 	}
    638     }
    639     return 1;
    640 }
    641 
    642 /*--------------------------------------------------------------------------
    643  * project - project a set of homogeneous coordinates into inhomogeneous ones
    644  *--------------------------------------------------------------------------
    645  */
    646 int
    647 Mapdesc::project( REAL *src, int stride, REAL *dest, int tstride, int ncols )
    648 {
    649     int s = sign( src[inhcoords] );
    650 
    651     REAL *clast = src + ncols * stride;
    652     for( REAL *cptr = src, *tcptr = dest; cptr != clast; cptr+=stride, tcptr+=tstride ) {
    653 	REAL *coordlast = cptr + inhcoords;
    654 	if( sign( *coordlast ) != s ) return 0;
    655 	for( REAL *coord = cptr, *tcoord = tcptr; coord != coordlast; coord++, tcoord++ )
    656 	    *tcoord = *coord / *coordlast;
    657     }
    658 
    659     return 1;
    660 }
    661 
    662 int
    663 Mapdesc::bboxTooBig(
    664     REAL *p,
    665     int	 rstride,
    666     int	 cstride,
    667     int	 nrows,
    668     int	 ncols,
    669     REAL bb[2][MAXCOORDS] )
    670 {
    671     REAL bbpts[MAXORDER][MAXORDER][MAXCOORDS];
    672     const int trstride = sizeof(bbpts[0]) / sizeof(REAL);
    673     const int tcstride = sizeof(bbpts[0][0]) / sizeof(REAL);
    674 
    675     // points have been transformed, therefore they are homogeneous
    676     // project points
    677     int val = project( p, rstride, cstride,
    678 	     &bbpts[0][0][0], trstride, tcstride, nrows, ncols );
    679     if( val == 0 ) return -1;
    680 
    681     // compute bounding box
    682     bbox( bb, &bbpts[0][0][0], trstride, tcstride, nrows, ncols );
    683 
    684     // find out if bounding box can't fit in unit cube
    685     if( bbox_subdividing == N_BBOXROUND ) {
    686 	for( int k=0; k != inhcoords; k++ )
    687 	    if( ceilf(bb[1][k]) - floorf(bb[0][k]) > bboxsize[k] ) return 1;
    688     } else {
    689 	for( int k=0; k != inhcoords; k++ )
    690 	    if( bb[1][k] - bb[0][k] > bboxsize[k] ) return 1;
    691     }
    692     return 0;
    693 }
    694 
    695 void
    696 Mapdesc::bbox(
    697     REAL bb[2][MAXCOORDS],
    698     REAL *p,
    699     int	 rstride,
    700     int	 cstride,
    701     int	 nrows,
    702     int	 ncols )
    703 {
    704     int k;
    705     for( k=0; k != inhcoords; k++ )
    706 	 bb[0][k] = bb[1][k] = p[k];
    707 
    708     for( int i=0; i != nrows; i++ )
    709 	for( int j=0; j != ncols; j++ )
    710 	    for( k=0; k != inhcoords; k++ ) {
    711 		REAL x = p[i*rstride + j*cstride + k];
    712 		if(  x < bb[0][k] ) bb[0][k] = x;
    713 		else if( x > bb[1][k] ) bb[1][k] = x;
    714 	    }
    715 }
    716 
    717 /*--------------------------------------------------------------------------
    718  * calcVelocityRational - calculate upper bound on first partial derivative
    719  *	of a homogeneous set of points and bounds on each row of points.
    720  *--------------------------------------------------------------------------
    721  */
    722 REAL
    723 Mapdesc::calcVelocityRational( REAL *p, int stride, int ncols )
    724 {
    725     REAL tmp[MAXORDER][MAXCOORDS];
    726 
    727     assert( ncols <= MAXORDER );
    728 
    729     const int tstride = sizeof(tmp[0]) / sizeof(REAL);
    730 
    731     if( project( p, stride, &tmp[0][0], tstride, ncols ) ) {
    732 	return calcPartialVelocity( &tmp[0][0], tstride, ncols, 1, 1.0 );
    733     } else { /* XXX */
    734 	return calcPartialVelocity( &tmp[0][0], tstride, ncols, 1, 1.0 );
    735     }
    736 }
    737 
    738 /*--------------------------------------------------------------------------
    739  * calcVelocityNonrational - calculate upper bound on  first partial
    740  *	derivative of a inhomogeneous set of points.
    741  *--------------------------------------------------------------------------
    742  */
    743 REAL
    744 Mapdesc::calcVelocityNonrational( REAL *pts, int stride, int ncols )
    745 {
    746     return calcPartialVelocity( pts, stride, ncols, 1, 1.0 );
    747 }
    748 
    749 int
    750 Mapdesc::isProperty( long property )
    751 {
    752     switch ( property ) {
    753 	case N_PIXEL_TOLERANCE:
    754 	case N_ERROR_TOLERANCE:
    755 	case N_CULLING:
    756 	case N_BBOX_SUBDIVIDING:
    757 	case N_S_STEPS:
    758 	case N_T_STEPS:
    759         case N_SAMPLINGMETHOD:
    760         case N_CLAMPFACTOR:
    761         case N_MINSAVINGS:
    762 	    return 1;
    763 	default:
    764 	    return 0;
    765     }
    766 }
    767 
    768 REAL
    769 Mapdesc::getProperty( long property )
    770 {
    771     switch ( property ) {
    772 	case N_PIXEL_TOLERANCE:
    773 	    return pixel_tolerance;
    774 	case N_ERROR_TOLERANCE:
    775 	    return error_tolerance;
    776 	case N_CULLING:
    777 	    return culling_method;
    778 	case N_BBOX_SUBDIVIDING:
    779 	    return bbox_subdividing;
    780 	case N_S_STEPS:
    781 	    return s_steps;
    782 	case N_T_STEPS:
    783 	    return t_steps;
    784         case N_SAMPLINGMETHOD:
    785 	    return sampling_method;
    786         case N_CLAMPFACTOR:
    787 	    return clampfactor;
    788         case N_MINSAVINGS:
    789 	    return minsavings;
    790 	default:
    791 	    abort();
    792 	    return -1; //not necessary, needed to shut up compiler
    793     }
    794 }
    795 
    796 void
    797 Mapdesc::setProperty( long property, REAL value )
    798 {
    799 
    800     switch ( property ) {
    801 	case N_PIXEL_TOLERANCE:
    802 	    pixel_tolerance = value;
    803 	    break;
    804 	case N_ERROR_TOLERANCE:
    805 	    error_tolerance = value;
    806 	    break;
    807 	case N_CULLING:
    808 	    culling_method = value;
    809 	    break;
    810 	case N_BBOX_SUBDIVIDING:
    811 	    if( value <= 0.0 ) value = N_NOBBOXSUBDIVISION;
    812 	    bbox_subdividing = value;
    813 	    break;
    814 	case N_S_STEPS:
    815 	    if( value < 0.0 ) value = 0.0;
    816 	    s_steps = value;
    817 	    maxrate = ( value < 0.0 ) ? 0.0 : value;
    818 	    maxsrate = ( value < 0.0 ) ? 0.0 : value;
    819 	    break;
    820 	case N_T_STEPS:
    821 	    if( value < 0.0 ) value = 0.0;
    822 	    t_steps = value;
    823 	    maxtrate = ( value < 0.0 ) ? 0.0 : value;
    824 	    break;
    825 	case N_SAMPLINGMETHOD:
    826 	    sampling_method = value;
    827 	    break;
    828 	case N_CLAMPFACTOR:
    829 	    if( value <= 0.0 ) value = N_NOCLAMPING;
    830 	    clampfactor = value;
    831 	    break;
    832 	case N_MINSAVINGS:
    833 	    if( value <= 0.0 ) value = N_NOSAVINGSSUBDIVISION;
    834 	    minsavings = value;
    835 	    break;
    836 	default:
    837 	    abort();
    838 	    break;
    839     }
    840 }
    841 
    842