1 /* ----------------------------------------------------------------------------
2 @COPYRIGHT :
3 Copyright 1993,1994,1995 David MacDonald,
4 McConnell Brain Imaging Centre,
5 Montreal Neurological Institute, McGill University.
6 Permission to use, copy, modify, and distribute this
7 software and its documentation for any purpose and without
8 fee is hereby granted, provided that the above copyright
9 notice appear in all copies. The author and McGill University
10 make no representations about the suitability of this
11 software for any purpose. It is provided "as is" without
12 express or implied warranty.
13 ---------------------------------------------------------------------------- */
14
15 #include <internal_volume_io.h>
16
17 /* ----------------------------- MNI Header -----------------------------------
18 @NAME : multiply_basis_matrices
19 @INPUT : n_derivs
20 n_degs
21 m1
22 m2
23 @OUTPUT : prod
24 @RETURNS :
25 @DESCRIPTION: Performs a matrix multiply of the basis matrix with the
26 powers of the u's positions. Steps through the
27 matrices in the appropriate strides. Could use the more
28 general multiply_matrices below, but is done this way for speed.
29 @METHOD :
30 @GLOBALS :
31 @CALLS :
32 @CREATED : Jan 21, 1995 David MacDonald
33 @MODIFIED :
34 ---------------------------------------------------------------------------- */
35
multiply_basis_matrices(int n_derivs,int n_degs,Real m1[],Real m2[],Real prod[])36 static void multiply_basis_matrices(
37 int n_derivs,
38 int n_degs,
39 Real m1[],
40 Real m2[],
41 Real prod[] )
42 {
43 int i, j, k, m2_inc;
44 Real *m1_ptr, *m1_ptr1;
45 Real *m2_ptr;
46 Real sum, *prod_ptr;
47
48 m1_ptr = m1;
49 prod_ptr = prod;
50 m2_inc = 1 - n_degs * n_degs;
51 for_less( i, 0, n_derivs )
52 {
53 m2_ptr = m2;
54 for_less( j, 0, n_degs )
55 {
56 sum = 0.0;
57 m1_ptr1 = m1_ptr;
58 for_less( k, 0, n_degs )
59 {
60 sum += (*m1_ptr1) * (*m2_ptr);
61 ++m1_ptr1;
62 m2_ptr += n_degs;
63 }
64 m2_ptr += m2_inc;
65 *prod_ptr = sum;
66 ++prod_ptr;
67 }
68
69 m1_ptr += n_degs;
70 }
71 }
72
73 /* ----------------------------- MNI Header -----------------------------------
74 @NAME : multiply_matrices
75 @INPUT : x1 - x size of first matrix
76 y1 - y size of first matrix
77 m1 - first matrix
78 sa1 - x stride of first matrix
79 sb1 - y stride of first matrix
80 : y2 - y size of second matrix (x size must be y1)
81 m2 - second matrix
82 sa2 - x stride of second matrix
83 sb2 - y stride of second matrix
84 sap - x stride of product matrix
85 sbp - y stride of product matrix
86 @OUTPUT : prod - product of m1 * m2
87 @RETURNS :
88 @DESCRIPTION: Multiplies the two matrices m1 and m2, placing the results in
89 prod.
90 @METHOD :
91 @GLOBALS :
92 @CALLS :
93 @CREATED : Jan. 21, 1995 David MacDonald
94 @MODIFIED :
95 ---------------------------------------------------------------------------- */
96
multiply_matrices(int x1,int y1,Real m1[],int sa1,int sb1,int y2,Real m2[],int sa2,int sb2,Real prod[],int sap,int sbp)97 static void multiply_matrices(
98 int x1,
99 int y1,
100 Real m1[],
101 int sa1,
102 int sb1,
103 int y2,
104 Real m2[],
105 int sa2,
106 int sb2,
107 Real prod[],
108 int sap,
109 int sbp )
110 {
111 int i, j, k;
112 Real *m1_ptr, *m1_ptr1, *m2_ptr;
113 Real sum, *prod_ptr;
114
115 m1_ptr = m1;
116 prod_ptr = prod;
117 sb2 -= y1 * sa2;
118 sap -= y2 * sbp;
119 for_less( i, 0, x1 )
120 {
121 m2_ptr = m2;
122 for_less( j, 0, y2 )
123 {
124 sum = 0.0;
125 m1_ptr1 = m1_ptr;
126 for_less( k, 0, y1 )
127 {
128 sum += (*m1_ptr1) * (*m2_ptr);
129 m1_ptr1 += sb1;
130 m2_ptr += sa2;
131 }
132
133 *prod_ptr = sum;
134
135 prod_ptr += sbp;
136 m2_ptr += sb2;
137 }
138 m1_ptr += sa1;
139 prod_ptr += sap;
140 }
141 }
142
143 #define MAX_DEGREE 4
144 #define MAX_DIMS 10
145 #define MAX_TOTAL_VALUES 4000
146
147 /* ----------------------------- MNI Header -----------------------------------
148 @NAME : spline_tensor_product
149 @INPUT : n_dims
150 positions[n_dims]
151 degrees[n_dims]
152 bases[n_dims][degrees[dim]*degrees[dim]]
153 n_values
154 coefs [n_values*degrees[0]*degrees[1]*...]
155 n_derivs[n_dims]
156 @OUTPUT : results[n_values*n_derivs[0]*n_derivs[1]*...]
157 @RETURNS :
158 @DESCRIPTION: Performs the spline tensor product necessary to evaluate.
159 Takes as input the number of dimensions, the position to
160 evaluate, the basis matrices defining the interpolation method,
161 and the control vertices (coefs). The resulting values and
162 derivatives are placed in the 1D array results, conceptually as a
163 (1+n_dims)-D array.
164 @METHOD :
165 @GLOBALS :
166 @CALLS :
167 @CREATED : Jan 21, 1995 David MacDonald
168 @MODIFIED :
169 ---------------------------------------------------------------------------- */
170
spline_tensor_product(int n_dims,Real positions[],int degrees[],Real * bases[],int n_values,Real coefs[],int n_derivs[],Real results[])171 VIOAPI void spline_tensor_product(
172 int n_dims,
173 Real positions[],
174 int degrees[],
175 Real *bases[],
176 int n_values,
177 Real coefs[],
178 int n_derivs[],
179 Real results[] )
180 {
181 int deriv, d, k, total_values, src;
182 int ind, prev_ind, max_degree, n_derivs_plus_1, deg;
183 int static_indices[MAX_DIMS];
184 int *indices, total_derivs;
185 Real *input_coefs, u_power, u;
186 Real static_us[MAX_DEGREE*MAX_DEGREE];
187 Real static_weights[MAX_DEGREE*MAX_DEGREE];
188 Real *us, *weights;
189 Real *tmp_results[2], *r;
190 Real static_tmp_results[2][MAX_TOTAL_VALUES];
191 BOOLEAN results_alloced;
192
193 /*--- check arguments */
194
195 max_degree = 2;
196 total_values = n_values;
197 total_derivs = 0;
198 for_less( d, 0, n_dims )
199 {
200 if( degrees[d] < 2 )
201 {
202 print_error(
203 "spline_tensor_product: Degree %d must be greater than 1.\n",
204 degrees[d] );
205 return;
206 }
207 if( degrees[d] > max_degree )
208 max_degree = degrees[d];
209 if( n_derivs[d] > total_derivs )
210 total_derivs = n_derivs[d];
211
212 total_values *= degrees[d];
213 }
214
215 /*--- determine if fixed size storage is large enough,
216 if not allocate memory */
217
218 if( n_dims > MAX_DIMS )
219 {
220 ALLOC( indices, n_dims );
221 }
222 else
223 {
224 indices = static_indices;
225 }
226
227 if( max_degree > MAX_DEGREE )
228 {
229 ALLOC( us, max_degree * max_degree );
230 ALLOC( weights, max_degree * max_degree );
231 }
232 else
233 {
234 us = static_us;
235 weights = static_weights;
236 }
237
238 if( total_values > MAX_TOTAL_VALUES )
239 {
240 ALLOC( tmp_results[0], total_values );
241 ALLOC( tmp_results[1], total_values );
242 results_alloced = TRUE;
243 }
244 else
245 {
246 tmp_results[0] = static_tmp_results[0];
247 tmp_results[1] = static_tmp_results[1];
248 results_alloced = FALSE;
249 }
250
251 input_coefs = coefs;
252
253 src = 0;
254
255 /*--- do each dimension */
256
257 for_less( d, 0, n_dims )
258 {
259 deg = degrees[d];
260 n_derivs_plus_1 = 1 + n_derivs[d];
261
262 /*--- fill in the top row of matrix of powers of u
263 = [1 u u^2 u^3 ...] for evaluating values */
264
265 u = positions[d];
266 u_power = 1.0;
267 us[0] = 1.0;
268 for_less( k, 1, deg )
269 {
270 u_power *= u;
271 us[k] = u_power;
272 }
273
274 /*--- fill in the rest of the n_derivs_plus_1 by degrees[d] matrix:
275 1 u u^2 u^3 ...
276 0 1 2u 3u^2 ...
277 0 0 2 6u ...
278 ... */
279
280 ind = deg;
281 for_less( deriv, 1, n_derivs_plus_1 )
282 {
283 for_less( k, 0, deriv )
284 {
285 us[ind] = 0.0;
286 ++ind;
287 }
288
289 prev_ind = IJ( deriv-1, deriv-1, deg );
290 for_less( k, deriv, deg )
291 {
292 us[ind] = us[prev_ind] * (Real) k;
293 ++ind;
294 ++prev_ind;
295 }
296 }
297
298 /*--- multiply the u's matrix by the spline basis to create weights */
299
300 multiply_basis_matrices( n_derivs_plus_1, deg, us, bases[d], weights );
301
302 total_values /= deg;
303
304 if( d == n_dims-1 )
305 r = results;
306 else
307 r = tmp_results[1-src];
308
309 /*--- multiply coefficient weights by the coefficients */
310
311 multiply_matrices( n_derivs_plus_1, deg, weights, deg, 1,
312 total_values, input_coefs, total_values, 1,
313 r, 1, n_derivs_plus_1 );
314
315 src = 1 - src;
316 input_coefs = tmp_results[src];
317
318 total_values *= n_derivs_plus_1;
319 }
320
321 /*--- check to free memory */
322
323 if( n_dims > MAX_DIMS )
324 {
325 FREE( indices );
326 }
327
328 if( max_degree > MAX_DEGREE )
329 {
330 FREE( us );
331 FREE( weights );
332 }
333
334 if( results_alloced )
335 {
336 FREE( tmp_results[0] );
337 FREE( tmp_results[1] );
338 }
339 }
340