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
calcPartialVelocity(REAL * p,int stride,int ncols,int partial,REAL range)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
calcPartialVelocity(REAL * dist,REAL * p,int rstride,int cstride,int nrows,int ncols,int spartial,int tpartial,REAL srange,REAL trange,int side)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