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 */ 53REAL 54Mapdesc::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 */ 104REAL 105Mapdesc::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