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