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