FreeWRL/FreeX3D  3.0.0
mapdesc.cc
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( dest ) );
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, register REAL alpha, register 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  register int nc = inhcoords;
208  register REAL pw = p[nc];
209  register REAL nw = -pw;
210  register 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