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
50Mapdesc::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
83void
84Mapdesc::setBboxsize( INREAL *mat )
85{
86    for( int i = 0; i != inhcoords; i++ )
87	bboxsize[i] = (REAL) mat[i];
88}
89
90void
91Mapdesc::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
98void
99Mapdesc::surfbbox( REAL bb[2][MAXCOORDS] )
100{
101    backend.surfbbox( type, bb[0], bb[1] );
102}
103
104void
105Mapdesc::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 */
118void
119Mapdesc::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 */
158void
159Mapdesc::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 */
201unsigned int
202Mapdesc::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 */
286void
287Mapdesc::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 */
320void
321Mapdesc::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 */
354int
355Mapdesc::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 */
392int
393Mapdesc::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 */
424int
425Mapdesc::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 */
453void
454Mapdesc::xformSampling( REAL *pts, int order, int stride, REAL *sp, int outstride )
455{
456    xformMat( smat, pts, order, stride, sp, outstride );
457}
458
459void
460Mapdesc::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 */
471void
472Mapdesc::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 */
483void
484Mapdesc::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 */
498void
499Mapdesc::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
507void
508Mapdesc::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
516void
517Mapdesc::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
540void
541Mapdesc::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
576void
577Mapdesc::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
594void
595Mapdesc::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 */
619int
620Mapdesc::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 */
646int
647Mapdesc::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
662int
663Mapdesc::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
695void
696Mapdesc::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 */
722REAL
723Mapdesc::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 */
743REAL
744Mapdesc::calcVelocityNonrational( REAL *pts, int stride, int ncols )
745{
746    return calcPartialVelocity( pts, stride, ncols, 1, 1.0 );
747}
748
749int
750Mapdesc::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
768REAL
769Mapdesc::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
796void
797Mapdesc::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