workbench/libs/mesa/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  */
  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[0][0] ) * 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, 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