1 /*
2 * vp_resample.c
3 *
4 * Routines to resample an array to a grid with a different resolution.
5 *
6 * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
7 * Junior University. All rights reserved.
8 *
9 * Permission to use, copy, modify and distribute this software and its
10 * documentation for any purpose is hereby granted without fee, provided
11 * that the above copyright notice and this permission notice appear in
12 * all copies of this software and that you do not sell the software.
13 * Commercial licensing is available by contacting the author.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
17 * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
18 *
19 * Author:
20 * Phil Lacroute
21 * Computer Systems Laboratory
22 * Electrical Engineering Dept.
23 * Stanford University
24 */
25
26 /*
27 * $Date: 1994/12/30 23:52:38 $
28 * $Revision: 1.6 $
29 */
30
31 #include "vp_global.h"
32
33 /* convert a float in the interval [0-1) to a 31-bit fixed point */
34 #define FLTFRAC_TO_FIX31(f) ((int)((f) * 2147483648.))
35
36 typedef struct {
37 int in_ptr_offset; /* offset in bytes from beginning of scanline
38 to first input sample for current
39 output sample */
40 float *wptr; /* filter weights for the filter phase
41 for current output sample */
42 int tap_min; /* first tap to evaluate */
43 int tap_max; /* last tap to evaluate */
44 } FilterTemplate;
45
46 static void ResampleUchar ANSI_ARGS((vpContext *vpc, int num_dimens,
47 int *src_dimens, int *dst_dimens, int *src_strides, int *dst_strides,
48 unsigned char *in_array, unsigned char *out_array,
49 FilterTemplate *template));
50 static void ResampleUshort ANSI_ARGS((vpContext *vpc, int num_dimens,
51 int *src_dimens, int *dst_dimens, int *src_strides, int *dst_strides,
52 unsigned short *in_array, unsigned short *out_array,
53 FilterTemplate *template));
54 static void ResampleFloat ANSI_ARGS((vpContext *vpc, int num_dimens,
55 int *src_dimens, int *dst_dimens, int *src_strides, int *dst_strides,
56 float *in_array, float *out_array, FilterTemplate *template));
57 static float *ComputeWeights ANSI_ARGS((vpContext *vpc, int src_xlen,
58 int dst_xlen, int filter_type));
59
60 /*
61 * vpSetFilter
62 *
63 * Set the filter to use for resampling.
64 */
65
66 vpResult
vpSetFilter(vpc,num_taps,num_phases,weights)67 vpSetFilter(vpc, num_taps, num_phases, weights)
68 vpContext *vpc;
69 int num_taps; /* number of filter taps */
70 int num_phases; /* number of filter phases */
71 float *weights; /* table of filter weights (weights[num_phases][num_taps]) */
72 {
73 int num_ones, bit;
74
75 /* make sure num_taps is positive and num_phases is a power of two */
76 if (num_taps < 1 || num_phases < 1)
77 return(VPSetError(vpc, VPERROR_BAD_VALUE));
78 num_ones = 0;
79 for (bit = 0; bit < 32; bit++) {
80 if (num_phases & (1 << bit))
81 num_ones++;
82 }
83 if (num_ones != 1)
84 return(VPSetError(vpc, VPERROR_BAD_VALUE));
85
86 /* store values in context */
87 vpc->filter_num_taps = num_taps;
88 vpc->filter_num_phases = num_phases;
89 vpc->filter_weights = weights;
90 return(VP_OK);
91 }
92
93 /*
94 * vpResample
95 *
96 * Resample an array to a grid with a different resolution.
97 */
98
99 vpResult
vpResample(vpc,num_dimens,src_dimens,dst_dimens,src_strides,dst_strides,element_type,in_array,out_array)100 vpResample(vpc, num_dimens, src_dimens, dst_dimens, src_strides, dst_strides,
101 element_type, in_array, out_array)
102 vpContext *vpc;
103 int num_dimens; /* number of dimensions in the two arrays */
104 int *src_dimens; /* sizes of source array dimensions */
105 int *dst_dimens; /* sizes of destination array dimensions (must
106 be the same, except for first dimension) */
107 int *src_strides; /* strides of source array dimensions (in bytes) */
108 int *dst_strides; /* strides of destination dimensions (in bytes) */
109 int element_type; /* type of array element (VP_UCHAR, VP_USHORT,
110 VP_FLOAT) */
111 void *in_array; /* input array containing data */
112 void *out_array; /* storage for output array */
113 {
114 int num_taps; /* number of filter taps */
115 int num_phases; /* number of filter phases */
116 int in_x_count; /* length of input scanlines */
117 int out_x_count; /* length of output scanlines */
118 int in_x_stride; /* stride of input scanline elements */
119 double scale_factor; /* in_x = scale_factor * out_x */
120 double in_x0; /* location of center of first output sample
121 in the input scanline */
122 int index0; /* coordinate of input sample corresponding
123 to first filter tap for first output
124 sample */
125 int phase0; /* filter phase for first output sample */
126 int index_incr; /* change in index0 for next output
127 sample */
128 int phase_incr; /* change in phase0 for next output
129 sample */
130 int unused_phase_bits; /* number of low-order bits of the phase that
131 are ignored for indexing the weight table */
132 FilterTemplate *template; /* filter template */
133 float *weights; /* pointer to weight table */
134 int in_offset; /* offset to input sample */
135 int index, phase; /* current input sample index and phase */
136 int out_x; /* current output sample */
137 int tap_min, tap_max; /* bounds on tap number */
138 int bit, d;
139
140 /* check for errors */
141 if (vpc->filter_weights == NULL)
142 return(VPSetError(vpc, VPERROR_BAD_SIZE));
143
144 /* find where the first output sample maps into the input array
145 and compute the filter phase for that sample; also compute
146 increments to get the input array position and filter phase
147 for the next sample */
148 num_taps = vpc->filter_num_taps;
149 num_phases = vpc->filter_num_phases;
150 in_x_count = src_dimens[0];
151 out_x_count = dst_dimens[0];
152 scale_factor = (double)in_x_count / (double)out_x_count;
153 if (num_taps % 2 == 0) {
154 /* even number of taps */
155
156 /* map center of first output voxel (x=0.5) to input voxel space
157 (multiply by scale_factor), then translate by -0.5 to convert
158 input voxels centered at 0.5 to input voxels centered at 0.0 */
159 in_x0 = 0.5 * scale_factor - 0.5;
160 phase0 = FLTFRAC_TO_FIX31(in_x0 - floor(in_x0));
161 index0 = (int)floor(in_x0) - num_taps/2 + 1;
162 } else {
163 /* odd number of taps */
164
165 /* omit translation by -0.5 since filter phase is offset by 0.5 voxels
166 relative to previous case */
167 in_x0 = 0.5 * scale_factor;
168 phase0 = FLTFRAC_TO_FIX31(in_x0 - floor(in_x0));
169 if (in_x0 < 0.5) {
170 index0 = (int)floor(in_x0) - num_taps/2;
171 } else {
172 index0 = (int)floor(in_x0) - num_taps/2 - 1;
173 }
174 }
175 index_incr = (int)floor(scale_factor);
176 phase_incr = FLTFRAC_TO_FIX31(scale_factor - index_incr);
177 unused_phase_bits = 0;
178 for (bit = 0; bit < 32; bit++) {
179 if (num_phases & (1 << bit)) {
180 unused_phase_bits = 31 - bit;
181 break;
182 }
183 }
184 ASSERT(unused_phase_bits != 0);
185
186 /* compute a template containing input array position and filter
187 weights for each output sample in an output scanline */
188 Alloc(vpc, template, FilterTemplate *, out_x_count*sizeof(FilterTemplate),
189 "FilterTemplate");
190 weights = vpc->filter_weights;
191 index = index0;
192 phase = phase0;
193 in_x_stride = src_strides[0];
194 in_offset = index * in_x_stride;
195 for (out_x = 0; out_x < out_x_count; out_x++) {
196 tap_min = MAX(0, -index);
197 tap_max = MIN(in_x_count - index - 1, num_taps-1);
198 template[out_x].in_ptr_offset = in_offset + tap_min * in_x_stride;
199 template[out_x].wptr = &weights[(phase >> unused_phase_bits) * num_taps
200 + tap_min];
201 template[out_x].tap_min = tap_min;
202 template[out_x].tap_max = tap_max;
203 phase += phase_incr;
204 if (phase < 0) {
205 phase &= 0x7FFFFFFF;
206 index += index_incr + 1;
207 in_offset += (index_incr + 1) * in_x_stride;
208 } else {
209 index += index_incr;
210 in_offset += index_incr * in_x_stride;
211 }
212 }
213
214 /* call a type-specific resampling routine */
215 switch (element_type) {
216 case VP_UCHAR:
217 ResampleUchar(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
218 dst_strides, in_array, out_array, template);
219 break;
220 case VP_USHORT:
221 ResampleUshort(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
222 dst_strides, in_array, out_array, template);
223 break;
224 case VP_FLOAT:
225 ResampleFloat(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
226 dst_strides, in_array, out_array, template);
227 break;
228 default:
229 Dealloc(vpc, template);
230 return(VPSetError(vpc, VPERROR_BAD_VALUE));
231 }
232 Dealloc(vpc, template);
233 return(VP_OK);
234 }
235
236 /*
237 * ResampleUchar
238 *
239 * Resample an array of unsigned chars.
240 */
241
242 static void
ResampleUchar(vpc,num_dimens,src_dimens,dst_dimens,src_strides,dst_strides,in_array,out_array,template)243 ResampleUchar(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
244 dst_strides, in_array, out_array, template)
245 vpContext *vpc;
246 int num_dimens; /* number of dimensions in the two arrays */
247 int *src_dimens; /* sizes of source array dimensions */
248 int *dst_dimens; /* sizes of destination array dimensions (must
249 be the same, except for first dimension) */
250 int *src_strides; /* strides of source array dimensions (in bytes) */
251 int *dst_strides; /* strides of destination dimensions (in bytes) */
252 unsigned char *in_array;/* input array containing data */
253 unsigned char *out_array;/* storage for output array */
254 FilterTemplate *template;/* filter template */
255 {
256 int out_x; /* current output sample */
257 float *wptr; /* pointer to filter weights */
258 float acc; /* accumulator for resampled value */
259 int tap; /* current tap number */
260 int tap_min, tap_max; /* bounds on tap number */
261 unsigned char *in_ptr; /* pointer to first input sample that
262 affects current output sample */
263 unsigned char *in_scan_ptr; /* pointer to beginning of input scanline */
264 unsigned char *out_ptr; /* pointer to current output sample */
265 unsigned char *out_scan_ptr;/* pointer to beginning of output scanline */
266 FilterTemplate *sample_template; /* template for output sample */
267 int out_x_count; /* number of elements in output scanline */
268 int in_x_stride; /* stride for input elements */
269 int out_x_stride; /* stride for output elements */
270 int *scan_coord; /* current scanline coordinates */
271 int done;
272 int dim;
273
274 /* copy parameters into local variables */
275 out_x_count = dst_dimens[0];
276 in_x_stride = src_strides[0];
277 out_x_stride = dst_strides[0];
278
279 /* allocate space for current scanline coordinates */
280 Alloc(vpc, scan_coord, int *, num_dimens * sizeof(int), "scan_coord");
281 for (dim = 0; dim < num_dimens; dim++) {
282 scan_coord[dim] = 0;
283 }
284
285 /* initialize pointers to first scanline */
286 in_scan_ptr = in_array;
287 out_scan_ptr = out_array;
288
289 done = 0;
290 while (!done) {
291 /* resample one scanline */
292 sample_template = template;
293 out_ptr = out_scan_ptr;
294 for (out_x = 0; out_x < out_x_count; out_x++) {
295 in_ptr = in_scan_ptr + sample_template->in_ptr_offset;
296 wptr = sample_template->wptr;
297 tap_min = sample_template->tap_min;
298 tap_max = sample_template->tap_max;
299 acc = 0;
300 for (tap = tap_min; tap <= tap_max; tap++) {
301 acc += (float)(*in_ptr) * *wptr;
302 in_ptr += in_x_stride;
303 wptr++;
304 }
305 if (acc > 255.)
306 *out_ptr = 255;
307 else if (acc < 0.)
308 *out_ptr = 0;
309 else
310 *out_ptr = (int)acc;
311 out_ptr += out_x_stride;
312 sample_template++;
313 } /* for out_x */
314
315 /* set pointers to next scanline */
316 for (dim = 1; dim < num_dimens; dim++) {
317 if (++scan_coord[dim] < src_dimens[dim]) {
318 in_scan_ptr += src_strides[dim];
319 out_scan_ptr += dst_strides[dim];
320 break;
321 } else if (dim == num_dimens-1) {
322 done = 1;
323 } else {
324 scan_coord[dim] = 0;
325 in_scan_ptr -= src_strides[dim] * src_dimens[dim];
326 out_scan_ptr -= dst_strides[dim] * dst_dimens[dim];
327 }
328 }
329 } /* while scanlines */
330
331 /* clean up */
332 Dealloc(vpc, scan_coord);
333 }
334
335 /*
336 * ResampleUshort
337 *
338 * Resample an array of unsigned shorts.
339 */
340
341 static void
ResampleUshort(vpc,num_dimens,src_dimens,dst_dimens,src_strides,dst_strides,in_array,out_array,template)342 ResampleUshort(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
343 dst_strides, in_array, out_array, template)
344 vpContext *vpc;
345 int num_dimens; /* number of dimensions in the two arrays */
346 int *src_dimens; /* sizes of source array dimensions */
347 int *dst_dimens; /* sizes of destination array dimensions (must
348 be the same, except for first dimension) */
349 int *src_strides; /* strides of source array dimensions (in bytes) */
350 int *dst_strides; /* strides of destination dimensions (in bytes) */
351 unsigned short *in_array;/* input array containing data */
352 unsigned short *out_array;/* storage for output array */
353 FilterTemplate *template;/* filter template */
354 {
355 int out_x; /* current output sample */
356 float *wptr; /* pointer to filter weights */
357 float acc; /* accumulator for resampled value */
358 int tap; /* current tap number */
359 int tap_min, tap_max; /* bounds on tap number */
360 unsigned short *in_ptr; /* pointer to first input sample that
361 affects current output sample */
362 unsigned short *in_scan_ptr;/* pointer to beginning of input scanline */
363 unsigned short *out_ptr; /* pointer to current output sample */
364 unsigned short *out_scan_ptr;/* pointer to beginning of output scanline */
365 FilterTemplate *sample_template; /* template for output sample */
366 int out_x_count; /* number of elements in output scanline */
367 int in_x_stride; /* stride for input elements */
368 int out_x_stride; /* stride for output elements */
369 int *scan_coord; /* current scanline coordinates */
370 int done;
371 int dim;
372
373 /* copy parameters into local variables */
374 out_x_count = dst_dimens[0];
375 in_x_stride = src_strides[0];
376 out_x_stride = dst_strides[0];
377
378 /* allocate space for current scanline coordinates */
379 Alloc(vpc, scan_coord, int *, num_dimens * sizeof(int), "scan_coord");
380 for (dim = 0; dim < num_dimens; dim++) {
381 scan_coord[dim] = 0;
382 }
383
384 /* initialize pointers to first scanline */
385 in_scan_ptr = in_array;
386 out_scan_ptr = out_array;
387
388 done = 0;
389 while (!done) {
390 /* resample one scanline */
391 sample_template = template;
392 out_ptr = out_scan_ptr;
393 for (out_x = 0; out_x < out_x_count; out_x++) {
394 in_ptr = in_scan_ptr + sample_template->in_ptr_offset;
395 wptr = sample_template->wptr;
396 tap_min = sample_template->tap_min;
397 tap_max = sample_template->tap_max;
398 acc = 0;
399 for (tap = tap_min; tap <= tap_max; tap++) {
400 acc += (float)(*in_ptr) * *wptr;
401 in_ptr = (unsigned short *)((char *)in_ptr + in_x_stride);
402 wptr++;
403 }
404 if (acc > 65535.)
405 *out_ptr = 65535;
406 else if (acc < 0.)
407 *out_ptr = 0;
408 else
409 *out_ptr = (int)acc;
410 out_ptr = (unsigned short *)((char *)out_ptr + out_x_stride);
411 sample_template++;
412 } /* for out_x */
413
414 /* set pointers to next scanline */
415 for (dim = 1; dim < num_dimens; dim++) {
416 if (++scan_coord[dim] < src_dimens[dim]) {
417 in_scan_ptr = (unsigned short *)((char *)in_scan_ptr +
418 src_strides[dim]);
419 out_scan_ptr = (unsigned short *)((char *)out_scan_ptr +
420 dst_strides[dim]);
421 break;
422 } else if (dim == num_dimens-1) {
423 done = 1;
424 } else {
425 scan_coord[dim] = 0;
426 in_scan_ptr = (unsigned short *)((char *)in_scan_ptr -
427 src_strides[dim] * src_dimens[dim]);
428 out_scan_ptr = (unsigned short *)((char *)out_scan_ptr -
429 dst_strides[dim] * dst_dimens[dim]);
430 }
431 }
432 } /* while scanlines */
433
434 /* clean up */
435 Dealloc(vpc, scan_coord);
436 }
437
438 /*
439 * ResampleFloat
440 *
441 * Resample an array of unsigned shorts.
442 */
443
444 static void
ResampleFloat(vpc,num_dimens,src_dimens,dst_dimens,src_strides,dst_strides,in_array,out_array,template)445 ResampleFloat(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
446 dst_strides, in_array, out_array, template)
447 vpContext *vpc;
448 int num_dimens; /* number of dimensions in the two arrays */
449 int *src_dimens; /* sizes of source array dimensions */
450 int *dst_dimens; /* sizes of destination array dimensions (must
451 be the same, except for first dimension) */
452 int *src_strides; /* strides of source array dimensions (in bytes) */
453 int *dst_strides; /* strides of destination dimensions (in bytes) */
454 float *in_array; /* input array containing data */
455 float *out_array; /* storage for output array */
456 FilterTemplate *template;/* filter template */
457 {
458 int out_x; /* current output sample */
459 float *wptr; /* pointer to filter weights */
460 float acc; /* accumulator for resampled value */
461 int tap; /* current tap number */
462 int tap_min, tap_max; /* bounds on tap number */
463 float *in_ptr; /* pointer to first input sample that
464 affects current output sample */
465 float *in_scan_ptr; /* pointer to beginning of input scanline */
466 float *out_ptr; /* pointer to current output sample */
467 float *out_scan_ptr; /* pointer to beginning of output scanline */
468 FilterTemplate *sample_template; /* template for output sample */
469 int out_x_count; /* number of elements in output scanline */
470 int in_x_stride; /* stride for input elements */
471 int out_x_stride; /* stride for output elements */
472 int *scan_coord; /* current scanline coordinates */
473 int done;
474 int dim;
475
476 /* copy parameters into local variables */
477 out_x_count = dst_dimens[0];
478 in_x_stride = src_strides[0];
479 out_x_stride = dst_strides[0];
480
481 /* allocate space for current scanline coordinates */
482 Alloc(vpc, scan_coord, int *, num_dimens * sizeof(int), "scan_coord");
483 for (dim = 0; dim < num_dimens; dim++) {
484 scan_coord[dim] = 0;
485 }
486
487 /* initialize pointers to first scanline */
488 in_scan_ptr = in_array;
489 out_scan_ptr = out_array;
490
491 done = 0;
492 while (!done) {
493 /* resample one scanline */
494 sample_template = template;
495 out_ptr = out_scan_ptr;
496 for (out_x = 0; out_x < out_x_count; out_x++) {
497 in_ptr = in_scan_ptr + sample_template->in_ptr_offset;
498 wptr = sample_template->wptr;
499 tap_min = sample_template->tap_min;
500 tap_max = sample_template->tap_max;
501 acc = 0;
502 for (tap = tap_min; tap <= tap_max; tap++) {
503 acc += *in_ptr * *wptr;
504 in_ptr = (float *)((char *)in_ptr + in_x_stride);
505 wptr++;
506 }
507 *out_ptr = acc;
508 out_ptr = (float *)((char *)out_ptr + out_x_stride);
509 sample_template++;
510 } /* for out_x */
511
512 /* set pointers to next scanline */
513 for (dim = 1; dim < num_dimens; dim++) {
514 if (++scan_coord[dim] < src_dimens[dim]) {
515 in_scan_ptr = (float *)((char *)in_scan_ptr +
516 src_strides[dim]);
517 out_scan_ptr = (float *)((char *)out_scan_ptr +
518 dst_strides[dim]);
519 break;
520 } else if (dim == num_dimens-1) {
521 done = 1;
522 } else {
523 scan_coord[dim] = 0;
524 in_scan_ptr = (float *)((char *)in_scan_ptr -
525 src_strides[dim] * src_dimens[dim]);
526 out_scan_ptr = (float *)((char *)out_scan_ptr -
527 dst_strides[dim] * dst_dimens[dim]);
528 }
529 }
530 } /* while scanlines */
531
532 /* clean up */
533 Dealloc(vpc, scan_coord);
534 }
535
536 /*
537 * vpResample2D
538 *
539 * Resample a 2D array.
540 */
541
542 vpResult
vpResample2D(in_array,in_x,in_y,out_array,out_x,out_y,element_type,filter_type)543 vpResample2D(in_array, in_x, in_y, out_array, out_x, out_y,
544 element_type, filter_type)
545 void *in_array; /* input array containing data */
546 int in_x, in_y; /* input array dimensions */
547 void *out_array; /* storage for output array */
548 int out_x, out_y; /* output array dimensions */
549 int element_type; /* type of array element (VP_UCHAR, VP_USHORT,
550 VP_FLOAT) */
551 int filter_type; /* type of filter (VP_BOX_FILTER, etc.) */
552 {
553 int src_dimens[2], dst_dimens[2];
554 int src_strides[2], dst_strides[2];
555 void *tmp1_array;
556 int element_size;
557 vpResult code;
558 vpContext *vpc;
559 float *weights;
560
561 /* compute size of array element and allocate intermediate arrays */
562 switch (element_type) {
563 case VP_UCHAR:
564 element_size = 1;
565 break;
566 case VP_USHORT:
567 element_size = 2;
568 break;
569 case VP_FLOAT:
570 element_size = 4;
571 break;
572 default:
573 return(VPSetError(vpc, VPERROR_BAD_OPTION));
574 }
575 vpc = vpCreateContext();
576 Alloc(vpc, tmp1_array, void *, out_x*in_y*element_size, "resample_tmp1");
577
578 /* resample first dimension */
579 src_dimens[0] = in_x;
580 src_dimens[1] = in_y;
581
582 dst_dimens[0] = out_x;
583 dst_dimens[1] = in_y;
584
585 src_strides[0] = element_size;
586 src_strides[1] = src_dimens[0] * src_strides[0];
587
588 dst_strides[0] = element_size;
589 dst_strides[1] = dst_dimens[0] * dst_strides[0];
590
591 weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
592 if (weights == NULL) {
593 Dealloc(vpc, tmp1_array);
594 return(vpc->error_code);
595 }
596 code = vpResample(vpc, 2, src_dimens, dst_dimens, src_strides, dst_strides,
597 element_type, in_array, tmp1_array);
598 Dealloc(vpc, weights);
599
600 if (code != VP_OK) {
601 Dealloc(vpc, tmp1_array);
602 return(code);
603 }
604
605 /* resample second dimension */
606 src_dimens[1] = out_x;
607 src_dimens[0] = in_y;
608
609 dst_dimens[1] = out_x;
610 dst_dimens[0] = out_y;
611
612 src_strides[1] = element_size;
613 src_strides[0] = src_dimens[1] * src_strides[1];
614
615 dst_strides[1] = element_size;
616 dst_strides[0] = dst_dimens[1] * dst_strides[1];
617
618 weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
619 if (weights == NULL) {
620 Dealloc(vpc, tmp1_array);
621 return(vpc->error_code);
622 }
623 code = vpResample(vpc, 2, src_dimens, dst_dimens, src_strides, dst_strides,
624 element_type, tmp1_array, out_array);
625 Dealloc(vpc, weights);
626
627 if (code != VP_OK) {
628 Dealloc(vpc, tmp1_array);
629 return(code);
630 }
631
632 /* clean up */
633 Dealloc(vpc, tmp1_array);
634 return(VP_OK);
635 }
636
637 /*
638 * vpResample3D
639 *
640 * Resample a 3D array.
641 */
642
643 vpResult
vpResample3D(in_array,in_x,in_y,in_z,out_array,out_x,out_y,out_z,element_type,filter_type)644 vpResample3D(in_array, in_x, in_y, in_z, out_array, out_x, out_y, out_z,
645 element_type, filter_type)
646 void *in_array; /* input array containing data */
647 int in_x, in_y, in_z; /* input array dimensions */
648 void *out_array; /* storage for output array */
649 int out_x, out_y, out_z;/* output array dimensions */
650 int element_type; /* type of array element (VP_UCHAR, VP_USHORT,
651 VP_FLOAT) */
652 int filter_type; /* type of filter (VP_BOX_FILTER, etc.) */
653 {
654 int src_dimens[3], dst_dimens[3];
655 int src_strides[3], dst_strides[3];
656 void *tmp1_array, *tmp2_array;
657 int element_size;
658 vpResult code;
659 vpContext *vpc;
660 float *weights;
661
662 /* compute size of array element and allocate intermediate arrays */
663 switch (element_type) {
664 case VP_UCHAR:
665 element_size = 1;
666 break;
667 case VP_USHORT:
668 element_size = 2;
669 break;
670 case VP_FLOAT:
671 element_size = 4;
672 break;
673 default:
674 return(VPSetError(vpc, VPERROR_BAD_OPTION));
675 }
676 vpc = vpCreateContext();
677 Alloc(vpc, tmp1_array, void *, out_x * in_y * in_z * element_size,
678 "resample_tmp1");
679 Alloc(vpc, tmp2_array, void *, out_x * out_y * in_z * element_size,
680 "resample_tmp2");
681
682 /* resample first dimension */
683 src_dimens[0] = in_x;
684 src_dimens[1] = in_y;
685 src_dimens[2] = in_z;
686
687 dst_dimens[0] = out_x;
688 dst_dimens[1] = in_y;
689 dst_dimens[2] = in_z;
690
691 src_strides[0] = element_size;
692 src_strides[1] = src_dimens[0] * src_strides[0];
693 src_strides[2] = src_dimens[1] * src_strides[1];
694
695 dst_strides[0] = element_size;
696 dst_strides[1] = dst_dimens[0] * dst_strides[0];
697 dst_strides[2] = dst_dimens[1] * dst_strides[1];
698
699 weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
700 if (weights == NULL) {
701 Dealloc(vpc, tmp1_array);
702 Dealloc(vpc, tmp2_array);
703 return(vpc->error_code);
704 }
705 code = vpResample(vpc, 3, src_dimens, dst_dimens, src_strides, dst_strides,
706 element_type, in_array, tmp1_array);
707 Dealloc(vpc, weights);
708
709 if (code != VP_OK) {
710 Dealloc(vpc, tmp1_array);
711 Dealloc(vpc, tmp2_array);
712 return(code);
713 }
714
715 /* resample second dimension */
716 src_dimens[1] = out_x;
717 src_dimens[0] = in_y;
718 src_dimens[2] = in_z;
719
720 dst_dimens[1] = out_x;
721 dst_dimens[0] = out_y;
722 dst_dimens[2] = in_z;
723
724 src_strides[1] = element_size;
725 src_strides[0] = src_dimens[1] * src_strides[1];
726 src_strides[2] = src_dimens[0] * src_strides[0];
727
728 dst_strides[1] = element_size;
729 dst_strides[0] = dst_dimens[1] * dst_strides[1];
730 dst_strides[2] = dst_dimens[0] * dst_strides[0];
731
732 weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
733 if (weights == NULL) {
734 Dealloc(vpc, tmp1_array);
735 Dealloc(vpc, tmp2_array);
736 return(vpc->error_code);
737 }
738 code = vpResample(vpc, 3, src_dimens, dst_dimens, src_strides, dst_strides,
739 element_type, tmp1_array, tmp2_array);
740 Dealloc(vpc, weights);
741
742 if (code != VP_OK) {
743 Dealloc(vpc, tmp1_array);
744 Dealloc(vpc, tmp2_array);
745 return(code);
746 }
747
748 /* resample third dimension */
749 src_dimens[1] = out_x;
750 src_dimens[2] = out_y;
751 src_dimens[0] = in_z;
752
753 dst_dimens[1] = out_x;
754 dst_dimens[2] = out_y;
755 dst_dimens[0] = out_z;
756
757 src_strides[1] = element_size;
758 src_strides[2] = src_dimens[1] * src_strides[1];
759 src_strides[0] = src_dimens[2] * src_strides[2];
760
761 dst_strides[1] = element_size;
762 dst_strides[2] = dst_dimens[1] * dst_strides[1];
763 dst_strides[0] = dst_dimens[2] * dst_strides[2];
764
765 weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
766 if (weights == NULL) {
767 Dealloc(vpc, tmp1_array);
768 Dealloc(vpc, tmp2_array);
769 return(vpc->error_code);
770 }
771 code = vpResample(vpc, 3, src_dimens, dst_dimens, src_strides, dst_strides,
772 element_type, tmp2_array, out_array);
773 Dealloc(vpc, weights);
774
775 if (code != VP_OK) {
776 Dealloc(vpc, tmp1_array);
777 Dealloc(vpc, tmp2_array);
778 return(code);
779 }
780
781 /* clean up */
782 Dealloc(vpc, tmp1_array);
783 Dealloc(vpc, tmp2_array);
784 return(VP_OK);
785 }
786
787 /*
788 * ComputeWeights
789 *
790 * Allocate and compute a filter weight table for a predefined filter type.
791 */
792
793 static float *
ComputeWeights(vpc,src_xlen,dst_xlen,filter_type)794 ComputeWeights(vpc, src_xlen, dst_xlen, filter_type)
795 vpContext *vpc; /* context for storing table */
796 int src_xlen; /* number of samples in input scanline */
797 int dst_xlen; /* number of samples in output scanline */
798 int filter_type; /* type of filter (VP_BOX_FILTER, etc.) */
799 {
800 double scale_factor;
801 int num_phases, num_taps, support, tap_limit, phases, table_size;
802 int code;
803 float *weights;
804
805 switch (filter_type) {
806 case VP_BOX_FILTER:
807 support = 1;
808 break;
809 case VP_LINEAR_FILTER:
810 support = 2;
811 break;
812 case VP_GAUSSIAN_FILTER:
813 support = 3;
814 break;
815 case VP_BSPLINE_FILTER:
816 case VP_MITCHELL_FILTER:
817 support = 4;
818 break;
819 default:
820 VPSetError(vpc, VPERROR_BAD_OPTION);
821 return(NULL);
822 }
823 scale_factor = (double)dst_xlen / (double)src_xlen;
824 if (scale_factor >= 1.0) {
825 num_taps = support;
826 num_phases = 1024;
827 } else {
828 num_taps = (double)support / scale_factor;
829 tap_limit = 4;
830 phases = 1024;
831 while (1) {
832 if (num_taps <= tap_limit) {
833 num_phases = phases;
834 break;
835 }
836 tap_limit *= 2;
837 phases /= 2;
838 if (phases <= 1) {
839 num_phases = 1;
840 break;
841 }
842 }
843 }
844 table_size = num_taps * num_phases * sizeof(float);
845 Alloc(vpc, weights, float *, table_size, "weight_table");
846 switch (filter_type) {
847 case VP_BOX_FILTER:
848 code = vpBoxFilter(num_taps, num_phases, weights, table_size);
849 if (code != VP_OK) {
850 Dealloc(vpc, weights);
851 VPSetError(vpc, code);
852 return(NULL);
853 }
854 break;
855 case VP_LINEAR_FILTER:
856 code = vpLinearFilter(num_taps, num_phases, weights, table_size);
857 if (code != VP_OK) {
858 Dealloc(vpc, weights);
859 VPSetError(vpc, code);
860 return(NULL);
861 }
862 break;
863 case VP_GAUSSIAN_FILTER:
864 code = vpGaussianFilter(VP_GAUSSIAN_SIGMA, num_taps, num_phases,
865 weights, table_size);
866 if (code != VP_OK) {
867 Dealloc(vpc, weights);
868 VPSetError(vpc, code);
869 return(NULL);
870 }
871 break;
872 case VP_BSPLINE_FILTER:
873 code = vpBicubicFilter(1.0, 0.0, num_taps, num_phases, weights,
874 table_size);
875 if (code != VP_OK) {
876 Dealloc(vpc, weights);
877 VPSetError(vpc, code);
878 return(NULL);
879 }
880 break;
881 case VP_MITCHELL_FILTER:
882 code = vpBicubicFilter(1.0/3.0, 1.0/3.0, num_taps, num_phases,
883 weights, table_size);
884 if (code != VP_OK) {
885 Dealloc(vpc, weights);
886 VPSetError(vpc, code);
887 return(NULL);
888 }
889 break;
890 }
891 vpSetFilter(vpc, num_taps, num_phases, weights);
892 return(weights);
893 }
894
895 /*
896 * vpBoxFilter
897 *
898 * Compute filter weights for box filter.
899 * For abs(x) < 0.5:
900 * k(x) = C
901 * (C is chosen so that k(x) integrates to 1).
902 * Otherwise:
903 * k(x) = 0
904 */
905
906 vpResult
vpBoxFilter(num_taps,num_phases,weights,weights_bytes)907 vpBoxFilter(num_taps, num_phases, weights, weights_bytes)
908 int num_taps; /* number of filter taps to compute */
909 int num_phases; /* number of phases to compute */
910 float *weights; /* array for storing filter weights
911 (num_taps*num_phases entries) */
912 int weights_bytes; /* size of array (for error checking) */
913 {
914 int p, t;
915 float *wptr;
916 double value;
917
918 if (weights_bytes != num_taps * num_phases * sizeof(float))
919 return(VPERROR_BAD_SIZE);
920 if (num_phases % 2 != 0)
921 return(VPERROR_BAD_VALUE);
922 wptr = weights;
923 value = 1.0 / (double)num_taps;
924 for (p = 0; p < num_phases; p++) {
925 for (t = 0; t < num_taps; t++) {
926 *wptr++ = value;
927 }
928 }
929 return(VP_OK);
930 }
931
932 /*
933 * vpLinearFilter
934 *
935 * Compute filter weights for linear interpolation.
936 * For abs(x) < 1:
937 * k(x) = C * (1 - abs(x))
938 * (C is chosen so that k(x) integrates to 1).
939 * Otherwise:
940 * k(x) = 0
941 */
942
943 vpResult
vpLinearFilter(num_taps,num_phases,weights,weights_bytes)944 vpLinearFilter(num_taps, num_phases, weights, weights_bytes)
945 int num_taps; /* number of filter taps to compute */
946 int num_phases; /* number of phases to compute */
947 float *weights; /* array for storing filter weights
948 (num_taps*num_phases entries) */
949 int weights_bytes; /* size of array (for error checking) */
950 {
951 int p, t;
952 float *wptr1, *wptr2;
953 double x0, delta_x, x, xa, tap_spacing, sum, normalize, value;
954
955 if (weights_bytes != num_taps * num_phases * sizeof(float))
956 return(VPERROR_BAD_SIZE);
957 if (num_phases % 2 != 0)
958 return(VPERROR_BAD_VALUE);
959 wptr1 = weights;
960 tap_spacing = 2.0 / (double)num_taps;
961 x0 = -tap_spacing * ((double)num_taps/2.0 - 1.0);
962 delta_x = tap_spacing / (double)num_phases;
963 for (p = 0; p < num_phases/2; p++) {
964 x = x0;
965 sum = 0;
966 for (t = 0; t < num_taps; t++) {
967 if (x < 0.0)
968 xa = -x;
969 else
970 xa = x;
971 value = 1.0 - xa;
972 wptr1[t] = value;
973 sum += value;
974 x += tap_spacing;
975 }
976 normalize = 1.0 / sum;
977 for (t = 0; t < num_taps; t++) {
978 wptr1[t] *= normalize;
979 }
980 wptr1 += num_taps;
981 x0 -= delta_x;
982 }
983 wptr2 = wptr1;
984 for (p = 0; p < num_phases/2; p++) {
985 for (t = 0; t < num_taps; t++) {
986 *wptr1++ = *--wptr2;
987 }
988 }
989 return(VP_OK);
990 }
991
992 /*
993 * vpBicubicFilter
994 *
995 * Compute filter weights for a Mitchell bicubic.
996 *
997 * See Mitchell, D.P. and Netravali, A.N., "Reconstruction filters in
998 * computer graphics," Proc. SIGGRAPH '88 (Computer Graphics V22 N4),
999 * p. 221-8.
1000 */
1001
1002 vpResult
vpBicubicFilter(b_value,c_value,num_taps,num_phases,weights,weights_bytes)1003 vpBicubicFilter(b_value, c_value, num_taps, num_phases, weights, weights_bytes)
1004 double b_value; /* b in the filter kernel equation */
1005 double c_value; /* c in the filter kernel equation */
1006 int num_taps; /* number of filter taps to compute */
1007 int num_phases; /* number of phases to compute */
1008 float *weights; /* array for storing filter weights
1009 (num_taps*num_phases entries) */
1010 int weights_bytes; /* size of array (for error checking) */
1011 {
1012 int p, t;
1013 float *wptr1, *wptr2;
1014 double x0, delta_x, x, xa, tap_spacing, sum, normalize, value;
1015
1016 if (weights_bytes != num_taps * num_phases * sizeof(float))
1017 return(VPERROR_BAD_SIZE);
1018 if (num_phases % 2 != 0)
1019 return(VPERROR_BAD_VALUE);
1020 wptr1 = weights;
1021 tap_spacing = 4.0 / (double)num_taps;
1022 x0 = -tap_spacing * ((double)num_taps/2.0 - 1.0);
1023 delta_x = tap_spacing / (double)num_phases;
1024 for (p = 0; p < num_phases/2; p++) {
1025 x = x0;
1026 sum = 0;
1027 for (t = 0; t < num_taps; t++) {
1028 if (x < 0.0)
1029 xa = -x;
1030 else
1031 xa = x;
1032 if (xa < 1.0) {
1033 value = (((12. - 9.*b_value - 6.*c_value)*xa - 18. +
1034 12.*b_value + 6.*c_value)*xa*xa + 6. -
1035 2.*b_value) * 1./6.;
1036 } else {
1037 value = ((((-b_value - 6.*c_value)*xa + 6.*b_value +
1038 30.*c_value)*xa - 12.*b_value - 48.*c_value)*xa +
1039 8.*b_value + 24.*c_value)* 1./6.;
1040 }
1041 wptr1[t] = value;
1042 sum += value;
1043 x += tap_spacing;
1044 }
1045 normalize = 1.0 / sum;
1046 for (t = 0; t < num_taps; t++) {
1047 wptr1[t] *= normalize;
1048 }
1049 wptr1 += num_taps;
1050 x0 -= delta_x;
1051 }
1052 wptr2 = wptr1;
1053 for (p = 0; p < num_phases/2; p++) {
1054 for (t = 0; t < num_taps; t++) {
1055 *wptr1++ = *--wptr2;
1056 }
1057 }
1058 return(VP_OK);
1059 }
1060
1061 /*
1062 * vpGaussianFilter
1063 *
1064 * Compute filter weights for a Gaussian.
1065 * For abs(x) <= 1.0:
1066 * k(x) = C * exp(-x*x/(2*sigma*sigma))
1067 * (C is chosen so that k(x) integrates to 1).
1068 * Otherwise:
1069 * k(x) = 0
1070 */
1071
1072 vpResult
vpGaussianFilter(sigma,num_taps,num_phases,weights,weights_bytes)1073 vpGaussianFilter(sigma, num_taps, num_phases, weights, weights_bytes)
1074 double sigma; /* standard deviation */
1075 int num_taps; /* number of filter taps to compute */
1076 int num_phases; /* number of phases to compute */
1077 float *weights; /* array for storing filter weights
1078 (num_taps*num_phases entries) */
1079 int weights_bytes; /* size of array (for error checking) */
1080 {
1081 int p, t;
1082 float *wptr1, *wptr2;
1083 double x0, delta_x, x, tap_spacing, sigma2_inv, sum, normalize, value;
1084
1085 if (weights_bytes != num_taps * num_phases * sizeof(float))
1086 return(VPERROR_BAD_SIZE);
1087 if (num_phases % 2 != 0)
1088 return(VPERROR_BAD_VALUE);
1089 wptr1 = weights;
1090 sigma2_inv = -1.0 / (2.0 * sigma * sigma);
1091 tap_spacing = 2.0 / (double)num_taps;
1092 x0 = -tap_spacing * ((double)num_taps/2.0 - 1.0);
1093 delta_x = tap_spacing / (double)num_phases;
1094 for (p = 0; p < num_phases/2; p++) {
1095 x = x0;
1096 sum = 0;
1097 for (t = 0; t < num_taps; t++) {
1098 value = exp(x*x*sigma2_inv);
1099 wptr1[t] = value;
1100 sum += value;
1101 x += tap_spacing;
1102 }
1103 normalize = 1.0 / sum;
1104 for (t = 0; t < num_taps; t++) {
1105 wptr1[t] *= normalize;
1106 }
1107 wptr1 += num_taps;
1108 x0 -= delta_x;
1109 }
1110 wptr2 = wptr1;
1111 for (p = 0; p < num_phases/2; p++) {
1112 for (t = 0; t < num_taps; t++) {
1113 *wptr1++ = *--wptr2;
1114 }
1115 }
1116 return(VP_OK);
1117 }
1118