src/glu/sgi/libnurbs/internals/mapdescv.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  * mapdescv.c++
  37  *
  38  * $Date: 2001/03/17 00:25:41 $ $Revision: 1.1 $
  39  * $Header: /home/krh/git/sync/mesa-cvs-repo/Mesa/src/glu/sgi/libnurbs/internals/mapdescv.cc,v 1.1 2001/03/17 00:25:41 brianp Exp $
  40  */
  41
  42 #include "glimports.h"
  43 #include "mystdio.h"
  44 #include "myassert.h"
  45 #include "mystring.h"
  46 #include "mymath.h"
  47 #include "nurbsconsts.h"
  48 #include "mapdesc.h"
  49
  50 /*--------------------------------------------------------------------------
  51  * calcPartialVelocity - calculate maximum magnitude of a given partial
  52  * derivative
  53  *--------------------------------------------------------------------------
  54  */
  55 REAL
  56 Mapdesc::calcPartialVelocity (
  57     REAL *p,
  58     int  stride,
  59     int  ncols,
  60     int  partial,
  61     REAL range )
  62 {
  63     REAL tmp[MAXORDER][MAXCOORDS];
  64     REAL mag[MAXORDER];
  65
  66     assert( ncols <= MAXORDER );
  67
  68     int j, k, t;
  69     // copy inhomogeneous control points into temporary array
  70     for( j=0; j != ncols; j++ )
  71         for( k=0; k != inhcoords; k++ )
  72             tmp[j][k] = p[j*stride + k];
  73
  74     for( t=0; t != partial; t++ )
  75         for( j=0; j != ncols-t-1; j++ )
  76             for( k=0; k != inhcoords; k++ )
  77                 tmp[j][k] = tmp[j+1][k] - tmp[j][k];
  78
  79     // compute magnitude and store in mag array
  80     for( j=0; j != ncols-partial; j++ ) {
  81         mag[j] = 0.0;
  82         for( k=0; k != inhcoords; k++ )
  83             mag[j] += tmp[j][k] * tmp[j][k];
  84     }
  85
  86     // compute scale factor
  87     REAL fac = 1;
  88     REAL invt = 1.0 / range;
  89     for( t = ncols-1; t != ncols-1-partial; t-- )
  90         fac *= t * invt;
  91
  92     // compute max magnitude of all entries in array
  93     REAL max = 0.0;
  94     for( j=0; j != ncols-partial; j++ )
  95         if( mag[j] > max ) max = mag[j];
  96     max = fac * ::sqrtf( (float) max );
  97
  98     return max;
  99 }
 100
 101 /*--------------------------------------------------------------------------
 102  * calcPartialVelocity - calculate maximum magnitude of a given partial
 103  * derivative
 104  *--------------------------------------------------------------------------
 105  */
 106 REAL
 107 Mapdesc::calcPartialVelocity (
 108     REAL *dist,
 109     REAL *p,
 110     int  rstride,
 111     int  cstride,
 112     int  nrows,
 113     int  ncols,
 114     int  spartial,
 115     int  tpartial,
 116     REAL srange,
 117     REAL trange,
 118     int  side )
 119 {
 120     REAL tmp[MAXORDER][MAXORDER][MAXCOORDS];
 121     REAL mag[MAXORDER][MAXORDER];
 122
 123     assert( nrows <= MAXORDER );
 124     assert( ncols <= MAXORDER );
 125
 126     REAL *tp = &tmp[0][0][0];
 127     REAL *mp = &mag[0][0];
 128     const int istride = sizeof( tmp[0]) / sizeof( tmp[0][0][0] );
 129     const int jstride = sizeof( tmp[0][0]) / sizeof( tmp[0][0][0] );
 130     const int kstride = sizeof( tmp[0][0][0]) / sizeof( tmp[0][0][0] );
 131     const int mistride = sizeof( mag[0]) / sizeof( mag[0][0] );
 132     const int mjstride = sizeof( mag[0][0]) / sizeof( mag[0][0] );
 133     const int idist = nrows * istride;
 134     const int jdist = ncols * jstride;
 135     const int kdist = inhcoords * kstride;
 136     const int id = idist - spartial * istride;
 137     const int jd = jdist - tpartial * jstride;
 138
 139     {
 140         // copy control points
 141         REAL *ti = tp;
 142         REAL *qi = p;
 143         REAL *til = tp + idist;
 144         for( ; ti != til; ) {
 145             REAL *tj = ti;
 146             REAL *qj = qi;
 147             REAL *tjl = ti + jdist;
 148             for( ; tj != tjl;  ) {
 149                 for( int k=0; k != inhcoords; k++ ) {
 150                     tj[k] = qj[k];
 151                 }
 152                 tj += jstride;
 153                 qj += cstride;
 154             }
 155             ti += istride;
 156             qi += rstride;
 157         }
 158     }
 159
 160     {
 161         // compute (s)-partial derivative control points
 162         REAL *til = tp + idist - istride;
 163         const REAL *till = til - ( spartial * istride );
 164         for( ; til != till; til -= istride )
 165             for( REAL *ti = tp; ti != til; ti += istride )
 166                 for( REAL *tj = ti, *tjl = tj + jdist; tj != tjl; tj += jstride )
 167                     for( int k=0; k != inhcoords; k++ )
 168                         tj[k] = tj[k+istride] - tj[k];
 169     }
 170
 171     {
 172         // compute (s,t)-partial derivative control points
 173         REAL *tjl = tp + jdist - jstride;
 174         const REAL *tjll = tjl - ( tpartial * jstride );
 175         for( ; tjl != tjll; tjl -= jstride )
 176             for( REAL *tj = tp; tj != tjl; tj += jstride )
 177                 for( REAL *ti = tj, *til = ti + id; ti != til; ti += istride )
 178                     for( int k=0; k != inhcoords; k++ )
 179                         ti[k] = ti[k+jstride] - ti[k];
 180
 181     }
 182
 183     REAL max = 0.0;
 184     {
 185         // compute magnitude and store in mag array
 186         memset( (void *) mp, 0, sizeof( mag ) );
 187         for( REAL *ti = tp, *mi = mp, *til = tp + id; ti != til; ti += istride, mi += mistride )
 188             for( REAL *tj = ti, *mj = mi, *tjl = ti + jd; tj != tjl; tj += jstride, mj += mjstride ) {
 189                 for( int k=0; k != inhcoords; k++ )
 190                    *mj += tj[k] * tj[k];
 191                 if( *mj > max ) max = *mj;
 192             }
 193
 194     }
 195
 196     int i, j;
 197
 198     // compute scale factor
 199     REAL fac = 1.0;
 200     {
 201         REAL invs = 1.0 / srange;
 202         REAL invt = 1.0 / trange;
 203         for( int s = nrows-1, slast = s-spartial; s != slast; s-- )
 204             fac *= s * invs;
 205         for( int t = ncols-1, tlast = t-tpartial; t != tlast; t-- )
 206             fac *= t * invt;
 207     }
 208
 209     if( side == 0 ) {
 210         // compute max magnitude of first and last column
 211         dist[0] = 0.0;
 212         dist[1] = 0.0;
 213         for( i=0; i != nrows-spartial; i++ ) {
 214             j = 0;
 215             if( mag[i][j] > dist[0] ) dist[0] = mag[i][j];
 216
 217             j = ncols-tpartial-1;
 218             if( mag[i][j] > dist[1] ) dist[1] = mag[i][j];
 219         }
 220         dist[0] = fac * ::sqrtf( dist[0] );
 221         dist[1] = fac * ::sqrtf( dist[1] );
 222     } else if( side == 1 ) {
 223         // compute max magnitude of first and last row
 224         dist[0] = 0.0;
 225         dist[1] = 0.0;
 226         for( j=0; j != ncols-tpartial; j++ ) {
 227             i = 0;
 228             if( mag[i][j] > dist[0] ) dist[0] = mag[i][j];
 229
 230             i = nrows-spartial-1;
 231             if( mag[i][j] > dist[1] ) dist[1] = mag[i][j];
 232         }
 233         dist[0] = fac * ::sqrtf( dist[0] );
 234         dist[1] = fac * ::sqrtf( dist[1] );
 235     }
 236
 237     max = fac * ::sqrtf( (float) max );
 238
 239     return max;
 240 }
 241