00001 /*
00002  * vp_resample.c
00003  *
00004  * Routines to resample an array to a grid with a different resolution.
00005  *
00006  * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
00007  * Junior University.  All rights reserved.
00008  *
00009  * Permission to use, copy, modify and distribute this software and its
00010  * documentation for any purpose is hereby granted without fee, provided
00011  * that the above copyright notice and this permission notice appear in
00012  * all copies of this software and that you do not sell the software.
00013  * Commercial licensing is available by contacting the author.
00014  * 
00015  * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
00016  * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
00017  * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
00018  *
00019  * Author:
00020  *    Phil Lacroute
00021  *    Computer Systems Laboratory
00022  *    Electrical Engineering Dept.
00023  *    Stanford University
00024  */
00025 
00026 /*
00027  * $Date: 2001/12/17 16:16:23 $
00028  * $Revision: 1.1 $
00029  */
00030 
00031 #include "vp_global.h"
00032 
00033 /* convert a float in the interval [0-1) to a 31-bit fixed point */
00034 #define FLTFRAC_TO_FIX31(f)     ((int)((f) * 2147483648.))
00035 
00036 typedef struct {
00037     int in_ptr_offset;          /* offset in bytes from beginning of scanline
00038                                    to first input sample for current 
00039                                    output sample */
00040     float *wptr;                /* filter weights for the filter phase
00041                                    for current output sample */
00042     int tap_min;                /* first tap to evaluate */
00043     int tap_max;                /* last tap to evaluate */
00044 } FilterTemplate;
00045 
00046 static void ResampleUchar ANSI_ARGS((vpContext *vpc, int num_dimens,
00047     int *src_dimens, int *dst_dimens, int *src_strides, int *dst_strides,
00048     unsigned char *in_array, unsigned char *out_array,
00049     FilterTemplate *template));
00050 static void ResampleUshort ANSI_ARGS((vpContext *vpc, int num_dimens,
00051     int *src_dimens, int *dst_dimens, int *src_strides, int *dst_strides,
00052     unsigned short *in_array, unsigned short *out_array,
00053     FilterTemplate *template));
00054 static void ResampleFloat ANSI_ARGS((vpContext *vpc, int num_dimens,
00055     int *src_dimens, int *dst_dimens, int *src_strides, int *dst_strides,
00056     float *in_array, float *out_array, FilterTemplate *template));
00057 static float *ComputeWeights ANSI_ARGS((vpContext *vpc, int src_xlen, 
00058     int dst_xlen, int filter_type));
00059 
00060 /*
00061  * vpSetFilter
00062  *
00063  * Set the filter to use for resampling.
00064  */
00065 
00066 vpResult
00067 vpSetFilter(vpc, num_taps, num_phases, weights)
00068 vpContext *vpc;
00069 int num_taps;   /* number of filter taps */
00070 int num_phases; /* number of filter phases */
00071 float *weights; /* table of filter weights (weights[num_phases][num_taps]) */
00072 {
00073     int num_ones, bit;
00074 
00075     /* make sure num_taps is positive and num_phases is a power of two */
00076     if (num_taps < 1 || num_phases < 1)
00077         return(VPSetError(vpc, VPERROR_BAD_VALUE));
00078     num_ones = 0;
00079     for (bit = 0; bit < 32; bit++) {
00080         if (num_phases & (1 << bit))
00081             num_ones++;
00082     }
00083     if (num_ones != 1)
00084         return(VPSetError(vpc, VPERROR_BAD_VALUE));
00085 
00086     /* store values in context */
00087     vpc->filter_num_taps = num_taps;
00088     vpc->filter_num_phases = num_phases;
00089     vpc->filter_weights = weights;
00090     return(VP_OK);
00091 }
00092 
00093 /*
00094  * vpResample
00095  *
00096  * Resample an array to a grid with a different resolution.
00097  */
00098 
00099 vpResult
00100 vpResample(vpc, num_dimens, src_dimens, dst_dimens, src_strides, dst_strides,
00101            element_type, in_array, out_array)
00102 vpContext *vpc;
00103 int num_dimens;         /* number of dimensions in the two arrays */
00104 int *src_dimens;        /* sizes of source array dimensions */
00105 int *dst_dimens;        /* sizes of destination array dimensions (must
00106                            be the same, except for first dimension) */
00107 int *src_strides;       /* strides of source array dimensions (in bytes) */
00108 int *dst_strides;       /* strides of destination dimensions (in bytes) */
00109 int element_type;       /* type of array element (VP_UCHAR, VP_USHORT,
00110                            VP_FLOAT) */
00111 void *in_array;         /* input array containing data */
00112 void *out_array;        /* storage for output array */
00113 {
00114     int num_taps;               /* number of filter taps */
00115     int num_phases;             /* number of filter phases */
00116     int in_x_count;             /* length of input scanlines */
00117     int out_x_count;            /* length of output scanlines */
00118     int in_x_stride;            /* stride of input scanline elements */
00119     double scale_factor;        /* in_x = scale_factor * out_x */
00120     double in_x0;               /* location of center of first output sample
00121                                    in the input scanline */
00122     int index0;                 /* coordinate of input sample corresponding
00123                                    to first filter tap for first output
00124                                    sample */
00125     int phase0;                 /* filter phase for first output sample */
00126     int index_incr;             /* change in index0 for next output
00127                                    sample */
00128     int phase_incr;             /* change in phase0 for next output
00129                                    sample */
00130     int unused_phase_bits;      /* number of low-order bits of the phase that
00131                                    are ignored for indexing the weight table */
00132     FilterTemplate *template;   /* filter template */
00133     float *weights;             /* pointer to weight table */
00134     int in_offset;              /* offset to input sample */
00135     int index, phase;           /* current input sample index and phase */
00136     int out_x;                  /* current output sample */
00137     int tap_min, tap_max;       /* bounds on tap number */
00138     int bit, d;
00139 
00140     /* check for errors */
00141     if (vpc->filter_weights == NULL)
00142         return(VPSetError(vpc, VPERROR_BAD_SIZE));
00143 
00144     /* find where the first output sample maps into the input array
00145        and compute the filter phase for that sample; also compute
00146        increments to get the input array position and filter phase
00147        for the next sample */
00148     num_taps = vpc->filter_num_taps;
00149     num_phases = vpc->filter_num_phases;
00150     in_x_count = src_dimens[0];
00151     out_x_count = dst_dimens[0];
00152     scale_factor = (double)in_x_count / (double)out_x_count;
00153     if (num_taps % 2 == 0) {
00154         /* even number of taps */
00155 
00156         /* map center of first output voxel (x=0.5) to input voxel space
00157            (multiply by scale_factor), then translate by -0.5 to convert
00158            input voxels centered at 0.5 to input voxels centered at 0.0 */
00159         in_x0 = 0.5 * scale_factor - 0.5;
00160         phase0 = FLTFRAC_TO_FIX31(in_x0 - floor(in_x0));
00161         index0 = (int)floor(in_x0) - num_taps/2 + 1;
00162     } else {
00163         /* odd number of taps */
00164 
00165         /* omit translation by -0.5 since filter phase is offset by 0.5 voxels
00166            relative to previous case */
00167         in_x0 = 0.5 * scale_factor;
00168         phase0 = FLTFRAC_TO_FIX31(in_x0 - floor(in_x0));
00169         if (in_x0 < 0.5) {
00170             index0 = (int)floor(in_x0) - num_taps/2;
00171         } else {
00172             index0 = (int)floor(in_x0) - num_taps/2 - 1;
00173         }
00174     }
00175     index_incr = (int)floor(scale_factor);
00176     phase_incr = FLTFRAC_TO_FIX31(scale_factor - index_incr);
00177     unused_phase_bits = 0;
00178     for (bit = 0; bit < 32; bit++) {
00179         if (num_phases & (1 << bit)) {
00180             unused_phase_bits = 31 - bit;
00181             break;
00182         }
00183     }
00184     ASSERT(unused_phase_bits != 0);
00185 
00186     /* compute a template containing input array position and filter
00187        weights for each output sample in an output scanline */
00188     Alloc(vpc, template, FilterTemplate *, out_x_count*sizeof(FilterTemplate),
00189           "FilterTemplate");
00190     weights = vpc->filter_weights;
00191     index = index0;
00192     phase = phase0;
00193     in_x_stride = src_strides[0];
00194     in_offset = index * in_x_stride;
00195     for (out_x = 0; out_x < out_x_count; out_x++) {
00196         tap_min = MAX(0, -index);
00197         tap_max = MIN(in_x_count - index - 1, num_taps-1);
00198         template[out_x].in_ptr_offset = in_offset + tap_min * in_x_stride;
00199         template[out_x].wptr = &weights[(phase >> unused_phase_bits) * num_taps
00200                                         + tap_min];
00201         template[out_x].tap_min = tap_min;
00202         template[out_x].tap_max = tap_max;
00203         phase += phase_incr;
00204         if (phase < 0) {
00205             phase &= 0x7FFFFFFF;
00206             index += index_incr + 1;
00207             in_offset += (index_incr + 1) * in_x_stride;
00208         } else {
00209             index += index_incr;
00210             in_offset += index_incr * in_x_stride;
00211         }
00212     }
00213 
00214     /* call a type-specific resampling routine */
00215     switch (element_type) {
00216     case VP_UCHAR:
00217         ResampleUchar(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
00218                       dst_strides, in_array, out_array, template);
00219         break;
00220     case VP_USHORT:
00221         ResampleUshort(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
00222                        dst_strides, in_array, out_array, template);
00223         break;
00224     case VP_FLOAT:
00225         ResampleFloat(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
00226                       dst_strides, in_array, out_array, template);
00227         break;
00228     default:
00229         Dealloc(vpc, template);
00230         return(VPSetError(vpc, VPERROR_BAD_VALUE));
00231     }
00232     Dealloc(vpc, template);
00233     return(VP_OK);
00234 }
00235 
00236 /*
00237  * ResampleUchar
00238  *
00239  * Resample an array of unsigned chars.
00240  */
00241 
00242 static void
00243 ResampleUchar(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
00244               dst_strides, in_array, out_array, template)
00245 vpContext *vpc;
00246 int num_dimens;         /* number of dimensions in the two arrays */
00247 int *src_dimens;        /* sizes of source array dimensions */
00248 int *dst_dimens;        /* sizes of destination array dimensions (must
00249                            be the same, except for first dimension) */
00250 int *src_strides;       /* strides of source array dimensions (in bytes) */
00251 int *dst_strides;       /* strides of destination dimensions (in bytes) */
00252 unsigned char *in_array;/* input array containing data */
00253 unsigned char *out_array;/* storage for output array */
00254 FilterTemplate *template;/* filter template */
00255 {
00256     int out_x;                  /* current output sample */
00257     float *wptr;                /* pointer to filter weights */
00258     float acc;                  /* accumulator for resampled value */
00259     int tap;                    /* current tap number */
00260     int tap_min, tap_max;       /* bounds on tap number */
00261     unsigned char *in_ptr;      /* pointer to first input sample that
00262                                    affects current output sample */
00263     unsigned char *in_scan_ptr; /* pointer to beginning of input scanline */
00264     unsigned char *out_ptr;     /* pointer to current output sample */
00265     unsigned char *out_scan_ptr;/* pointer to beginning of output scanline */
00266     FilterTemplate *sample_template;    /* template for output sample */
00267     int out_x_count;            /* number of elements in output scanline */
00268     int in_x_stride;            /* stride for input elements */
00269     int out_x_stride;           /* stride for output elements */
00270     int *scan_coord;            /* current scanline coordinates */
00271     int done;
00272     int dim;
00273 
00274     /* copy parameters into local variables */
00275     out_x_count = dst_dimens[0];
00276     in_x_stride = src_strides[0];
00277     out_x_stride = dst_strides[0];
00278     
00279     /* allocate space for current scanline coordinates */
00280     Alloc(vpc, scan_coord, int *, num_dimens * sizeof(int), "scan_coord");
00281     for (dim = 0; dim < num_dimens; dim++) {
00282         scan_coord[dim] = 0;
00283     }
00284 
00285     /* initialize pointers to first scanline */
00286     in_scan_ptr = in_array;
00287     out_scan_ptr = out_array;
00288 
00289     done = 0;
00290     while (!done) {
00291         /* resample one scanline */
00292         sample_template = template;
00293         out_ptr = out_scan_ptr;
00294         for (out_x = 0; out_x < out_x_count; out_x++) {
00295             in_ptr = in_scan_ptr + sample_template->in_ptr_offset;
00296             wptr = sample_template->wptr;
00297             tap_min = sample_template->tap_min;
00298             tap_max = sample_template->tap_max;
00299             acc = 0;
00300             for (tap = tap_min; tap <= tap_max; tap++) {
00301                 acc += (float)(*in_ptr) * *wptr;
00302                 in_ptr += in_x_stride;
00303                 wptr++;
00304             }
00305             if (acc > 255.)
00306                 *out_ptr = 255;
00307             else if (acc < 0.)
00308                 *out_ptr = 0;
00309             else
00310                 *out_ptr = (int)acc;
00311             out_ptr += out_x_stride;
00312             sample_template++;
00313         } /* for out_x */
00314 
00315         /* set pointers to next scanline */
00316         for (dim = 1; dim < num_dimens; dim++) {
00317             if (++scan_coord[dim] < src_dimens[dim]) {
00318                 in_scan_ptr += src_strides[dim];
00319                 out_scan_ptr += dst_strides[dim];
00320                 break;
00321             } else if (dim == num_dimens-1) {
00322                 done = 1;
00323             } else {
00324                 scan_coord[dim] = 0;
00325                 in_scan_ptr -= src_strides[dim] * src_dimens[dim];
00326                 out_scan_ptr -= dst_strides[dim] * dst_dimens[dim];
00327             }
00328         }
00329     } /* while scanlines */
00330 
00331     /* clean up */
00332     Dealloc(vpc, scan_coord);
00333 }
00334 
00335 /*
00336  * ResampleUshort
00337  *
00338  * Resample an array of unsigned shorts.
00339  */
00340 
00341 static void
00342 ResampleUshort(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
00343                dst_strides, in_array, out_array, template)
00344 vpContext *vpc;
00345 int num_dimens;         /* number of dimensions in the two arrays */
00346 int *src_dimens;        /* sizes of source array dimensions */
00347 int *dst_dimens;        /* sizes of destination array dimensions (must
00348                            be the same, except for first dimension) */
00349 int *src_strides;       /* strides of source array dimensions (in bytes) */
00350 int *dst_strides;       /* strides of destination dimensions (in bytes) */
00351 unsigned short *in_array;/* input array containing data */
00352 unsigned short *out_array;/* storage for output array */
00353 FilterTemplate *template;/* filter template */
00354 {
00355     int out_x;                  /* current output sample */
00356     float *wptr;                /* pointer to filter weights */
00357     float acc;                  /* accumulator for resampled value */
00358     int tap;                    /* current tap number */
00359     int tap_min, tap_max;       /* bounds on tap number */
00360     unsigned short *in_ptr;     /* pointer to first input sample that
00361                                    affects current output sample */
00362     unsigned short *in_scan_ptr;/* pointer to beginning of input scanline */
00363     unsigned short *out_ptr;    /* pointer to current output sample */
00364     unsigned short *out_scan_ptr;/* pointer to beginning of output scanline */
00365     FilterTemplate *sample_template;    /* template for output sample */
00366     int out_x_count;            /* number of elements in output scanline */
00367     int in_x_stride;            /* stride for input elements */
00368     int out_x_stride;           /* stride for output elements */
00369     int *scan_coord;            /* current scanline coordinates */
00370     int done;
00371     int dim;
00372 
00373     /* copy parameters into local variables */
00374     out_x_count = dst_dimens[0];
00375     in_x_stride = src_strides[0];
00376     out_x_stride = dst_strides[0];
00377     
00378     /* allocate space for current scanline coordinates */
00379     Alloc(vpc, scan_coord, int *, num_dimens * sizeof(int), "scan_coord");
00380     for (dim = 0; dim < num_dimens; dim++) {
00381         scan_coord[dim] = 0;
00382     }
00383 
00384     /* initialize pointers to first scanline */
00385     in_scan_ptr = in_array;
00386     out_scan_ptr = out_array;
00387 
00388     done = 0;
00389     while (!done) {
00390         /* resample one scanline */
00391         sample_template = template;
00392         out_ptr = out_scan_ptr;
00393         for (out_x = 0; out_x < out_x_count; out_x++) {
00394             in_ptr = in_scan_ptr + sample_template->in_ptr_offset;
00395             wptr = sample_template->wptr;
00396             tap_min = sample_template->tap_min;
00397             tap_max = sample_template->tap_max;
00398             acc = 0;
00399             for (tap = tap_min; tap <= tap_max; tap++) {
00400                 acc += (float)(*in_ptr) * *wptr;
00401                 in_ptr = (unsigned short *)((char *)in_ptr + in_x_stride);
00402                 wptr++;
00403             }
00404             if (acc > 65535.)
00405                 *out_ptr = 65535;
00406             else if (acc < 0.)
00407                 *out_ptr = 0;
00408             else
00409                 *out_ptr = (int)acc;
00410             out_ptr = (unsigned short *)((char *)out_ptr + out_x_stride);
00411             sample_template++;
00412         } /* for out_x */
00413 
00414         /* set pointers to next scanline */
00415         for (dim = 1; dim < num_dimens; dim++) {
00416             if (++scan_coord[dim] < src_dimens[dim]) {
00417                 in_scan_ptr = (unsigned short *)((char *)in_scan_ptr +
00418                                                  src_strides[dim]);
00419                 out_scan_ptr = (unsigned short *)((char *)out_scan_ptr +
00420                                                   dst_strides[dim]);
00421                 break;
00422             } else if (dim == num_dimens-1) {
00423                 done = 1;
00424             } else {
00425                 scan_coord[dim] = 0;
00426                 in_scan_ptr = (unsigned short *)((char *)in_scan_ptr -
00427                               src_strides[dim] * src_dimens[dim]);
00428                 out_scan_ptr = (unsigned short *)((char *)out_scan_ptr -
00429                               dst_strides[dim] * dst_dimens[dim]);
00430             }
00431         }
00432     } /* while scanlines */
00433 
00434     /* clean up */
00435     Dealloc(vpc, scan_coord);
00436 }
00437 
00438 /*
00439  * ResampleFloat
00440  *
00441  * Resample an array of unsigned shorts.
00442  */
00443 
00444 static void
00445 ResampleFloat(vpc, num_dimens, src_dimens, dst_dimens, src_strides,
00446               dst_strides, in_array, out_array, template)
00447 vpContext *vpc;
00448 int num_dimens;         /* number of dimensions in the two arrays */
00449 int *src_dimens;        /* sizes of source array dimensions */
00450 int *dst_dimens;        /* sizes of destination array dimensions (must
00451                            be the same, except for first dimension) */
00452 int *src_strides;       /* strides of source array dimensions (in bytes) */
00453 int *dst_strides;       /* strides of destination dimensions (in bytes) */
00454 float *in_array;        /* input array containing data */
00455 float *out_array;       /* storage for output array */
00456 FilterTemplate *template;/* filter template */
00457 {
00458     int out_x;                  /* current output sample */
00459     float *wptr;                /* pointer to filter weights */
00460     float acc;                  /* accumulator for resampled value */
00461     int tap;                    /* current tap number */
00462     int tap_min, tap_max;       /* bounds on tap number */
00463     float *in_ptr;              /* pointer to first input sample that
00464                                    affects current output sample */
00465     float *in_scan_ptr;         /* pointer to beginning of input scanline */
00466     float *out_ptr;             /* pointer to current output sample */
00467     float *out_scan_ptr;        /* pointer to beginning of output scanline */
00468     FilterTemplate *sample_template;    /* template for output sample */
00469     int out_x_count;            /* number of elements in output scanline */
00470     int in_x_stride;            /* stride for input elements */
00471     int out_x_stride;           /* stride for output elements */
00472     int *scan_coord;            /* current scanline coordinates */
00473     int done;
00474     int dim;
00475 
00476     /* copy parameters into local variables */
00477     out_x_count = dst_dimens[0];
00478     in_x_stride = src_strides[0];
00479     out_x_stride = dst_strides[0];
00480     
00481     /* allocate space for current scanline coordinates */
00482     Alloc(vpc, scan_coord, int *, num_dimens * sizeof(int), "scan_coord");
00483     for (dim = 0; dim < num_dimens; dim++) {
00484         scan_coord[dim] = 0;
00485     }
00486 
00487     /* initialize pointers to first scanline */
00488     in_scan_ptr = in_array;
00489     out_scan_ptr = out_array;
00490 
00491     done = 0;
00492     while (!done) {
00493         /* resample one scanline */
00494         sample_template = template;
00495         out_ptr = out_scan_ptr;
00496         for (out_x = 0; out_x < out_x_count; out_x++) {
00497             in_ptr = in_scan_ptr + sample_template->in_ptr_offset;
00498             wptr = sample_template->wptr;
00499             tap_min = sample_template->tap_min;
00500             tap_max = sample_template->tap_max;
00501             acc = 0;
00502             for (tap = tap_min; tap <= tap_max; tap++) {
00503                 acc += *in_ptr * *wptr;
00504                 in_ptr = (float *)((char *)in_ptr + in_x_stride);
00505                 wptr++;
00506             }
00507             *out_ptr = acc;
00508             out_ptr = (float *)((char *)out_ptr + out_x_stride);
00509             sample_template++;
00510         } /* for out_x */
00511 
00512         /* set pointers to next scanline */
00513         for (dim = 1; dim < num_dimens; dim++) {
00514             if (++scan_coord[dim] < src_dimens[dim]) {
00515                 in_scan_ptr = (float *)((char *)in_scan_ptr +
00516                                         src_strides[dim]);
00517                 out_scan_ptr = (float *)((char *)out_scan_ptr +
00518                                          dst_strides[dim]);
00519                 break;
00520             } else if (dim == num_dimens-1) {
00521                 done = 1;
00522             } else {
00523                 scan_coord[dim] = 0;
00524                 in_scan_ptr = (float *)((char *)in_scan_ptr -
00525                               src_strides[dim] * src_dimens[dim]);
00526                 out_scan_ptr = (float *)((char *)out_scan_ptr -
00527                               dst_strides[dim] * dst_dimens[dim]);
00528             }
00529         }
00530     } /* while scanlines */
00531 
00532     /* clean up */
00533     Dealloc(vpc, scan_coord);
00534 }
00535 
00536 /*
00537  * vpResample2D
00538  *
00539  * Resample a 2D array.
00540  */
00541 
00542 vpResult
00543 vpResample2D(in_array, in_x, in_y, out_array, out_x, out_y,
00544              element_type, filter_type)
00545 void *in_array;         /* input array containing data */
00546 int in_x, in_y;         /* input array dimensions */
00547 void *out_array;        /* storage for output array */
00548 int out_x, out_y;       /* output array dimensions */
00549 int element_type;       /* type of array element (VP_UCHAR, VP_USHORT,
00550                            VP_FLOAT) */
00551 int filter_type;        /* type of filter (VP_BOX_FILTER, etc.) */
00552 {
00553     int src_dimens[2], dst_dimens[2];
00554     int src_strides[2], dst_strides[2];
00555     void *tmp1_array;
00556     int element_size;
00557     vpResult code;
00558     vpContext *vpc;
00559     float *weights;
00560 
00561     /* compute size of array element and allocate intermediate arrays */
00562     switch (element_type) {
00563     case VP_UCHAR:
00564         element_size = 1;
00565         break;
00566     case VP_USHORT:
00567         element_size = 2;
00568         break;
00569     case VP_FLOAT:
00570         element_size = 4;
00571         break;
00572     default:
00573         return(VPSetError(vpc, VPERROR_BAD_OPTION));
00574     }
00575     vpc = vpCreateContext();
00576     Alloc(vpc, tmp1_array, void *, out_x*in_y*element_size, "resample_tmp1");
00577 
00578     /* resample first dimension */
00579     src_dimens[0] = in_x;
00580     src_dimens[1] = in_y;
00581 
00582     dst_dimens[0] = out_x;
00583     dst_dimens[1] = in_y;
00584 
00585     src_strides[0] = element_size;
00586     src_strides[1] = src_dimens[0] * src_strides[0];
00587 
00588     dst_strides[0] = element_size;
00589     dst_strides[1] = dst_dimens[0] * dst_strides[0];
00590 
00591     weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
00592     if (weights == NULL) {
00593         Dealloc(vpc, tmp1_array);
00594         return(vpc->error_code);
00595     }
00596     code = vpResample(vpc, 2, src_dimens, dst_dimens, src_strides, dst_strides,
00597                       element_type, in_array, tmp1_array);
00598     Dealloc(vpc, weights);
00599 
00600     if (code != VP_OK) {
00601         Dealloc(vpc, tmp1_array);
00602         return(code);
00603     }
00604 
00605     /* resample second dimension */
00606     src_dimens[1] = out_x;
00607     src_dimens[0] = in_y;
00608 
00609     dst_dimens[1] = out_x;
00610     dst_dimens[0] = out_y;
00611 
00612     src_strides[1] = element_size;
00613     src_strides[0] = src_dimens[1] * src_strides[1];
00614 
00615     dst_strides[1] = element_size;
00616     dst_strides[0] = dst_dimens[1] * dst_strides[1];
00617 
00618     weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
00619     if (weights == NULL) {
00620         Dealloc(vpc, tmp1_array);
00621         return(vpc->error_code);
00622     }
00623     code = vpResample(vpc, 2, src_dimens, dst_dimens, src_strides, dst_strides,
00624                       element_type, tmp1_array, out_array);
00625     Dealloc(vpc, weights);
00626 
00627     if (code != VP_OK) {
00628         Dealloc(vpc, tmp1_array);
00629         return(code);
00630     }
00631 
00632     /* clean up */
00633     Dealloc(vpc, tmp1_array);
00634     return(VP_OK);
00635 }
00636 
00637 /*
00638  * vpResample3D
00639  *
00640  * Resample a 3D array.
00641  */
00642 
00643 vpResult
00644 vpResample3D(in_array, in_x, in_y, in_z, out_array, out_x, out_y, out_z,
00645              element_type, filter_type)
00646 void *in_array;         /* input array containing data */
00647 int in_x, in_y, in_z;   /* input array dimensions */
00648 void *out_array;        /* storage for output array */
00649 int out_x, out_y, out_z;/* output array dimensions */
00650 int element_type;       /* type of array element (VP_UCHAR, VP_USHORT,
00651                            VP_FLOAT) */
00652 int filter_type;        /* type of filter (VP_BOX_FILTER, etc.) */
00653 {
00654     int src_dimens[3], dst_dimens[3];
00655     int src_strides[3], dst_strides[3];
00656     void *tmp1_array, *tmp2_array;
00657     int element_size;
00658     vpResult code;
00659     vpContext *vpc;
00660     float *weights;
00661 
00662     /* compute size of array element and allocate intermediate arrays */
00663     switch (element_type) {
00664     case VP_UCHAR:
00665         element_size = 1;
00666         break;
00667     case VP_USHORT:
00668         element_size = 2;
00669         break;
00670     case VP_FLOAT:
00671         element_size = 4;
00672         break;
00673     default:
00674         return(VPSetError(vpc, VPERROR_BAD_OPTION));
00675     }
00676     vpc = vpCreateContext();
00677     Alloc(vpc, tmp1_array, void *, out_x * in_y * in_z * element_size,
00678           "resample_tmp1");
00679     Alloc(vpc, tmp2_array, void *, out_x * out_y * in_z * element_size,
00680           "resample_tmp2");
00681 
00682     /* resample first dimension */
00683     src_dimens[0] = in_x;
00684     src_dimens[1] = in_y;
00685     src_dimens[2] = in_z;
00686 
00687     dst_dimens[0] = out_x;
00688     dst_dimens[1] = in_y;
00689     dst_dimens[2] = in_z;
00690 
00691     src_strides[0] = element_size;
00692     src_strides[1] = src_dimens[0] * src_strides[0];
00693     src_strides[2] = src_dimens[1] * src_strides[1];
00694 
00695     dst_strides[0] = element_size;
00696     dst_strides[1] = dst_dimens[0] * dst_strides[0];
00697     dst_strides[2] = dst_dimens[1] * dst_strides[1];
00698 
00699     weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
00700     if (weights == NULL) {
00701         Dealloc(vpc, tmp1_array);
00702         Dealloc(vpc, tmp2_array);
00703         return(vpc->error_code);
00704     }
00705     code = vpResample(vpc, 3, src_dimens, dst_dimens, src_strides, dst_strides,
00706                       element_type, in_array, tmp1_array);
00707     Dealloc(vpc, weights);
00708 
00709     if (code != VP_OK) {
00710         Dealloc(vpc, tmp1_array);
00711         Dealloc(vpc, tmp2_array);
00712         return(code);
00713     }
00714 
00715     /* resample second dimension */
00716     src_dimens[1] = out_x;
00717     src_dimens[0] = in_y;
00718     src_dimens[2] = in_z;
00719 
00720     dst_dimens[1] = out_x;
00721     dst_dimens[0] = out_y;
00722     dst_dimens[2] = in_z;
00723 
00724     src_strides[1] = element_size;
00725     src_strides[0] = src_dimens[1] * src_strides[1];
00726     src_strides[2] = src_dimens[0] * src_strides[0];
00727 
00728     dst_strides[1] = element_size;
00729     dst_strides[0] = dst_dimens[1] * dst_strides[1];
00730     dst_strides[2] = dst_dimens[0] * dst_strides[0];
00731 
00732     weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
00733     if (weights == NULL) {
00734         Dealloc(vpc, tmp1_array);
00735         Dealloc(vpc, tmp2_array);
00736         return(vpc->error_code);
00737     }
00738     code = vpResample(vpc, 3, src_dimens, dst_dimens, src_strides, dst_strides,
00739                       element_type, tmp1_array, tmp2_array);
00740     Dealloc(vpc, weights);
00741 
00742     if (code != VP_OK) {
00743         Dealloc(vpc, tmp1_array);
00744         Dealloc(vpc, tmp2_array);
00745         return(code);
00746     }
00747 
00748     /* resample third dimension */
00749     src_dimens[1] = out_x;
00750     src_dimens[2] = out_y;
00751     src_dimens[0] = in_z;
00752 
00753     dst_dimens[1] = out_x;
00754     dst_dimens[2] = out_y;
00755     dst_dimens[0] = out_z;
00756 
00757     src_strides[1] = element_size;
00758     src_strides[2] = src_dimens[1] * src_strides[1];
00759     src_strides[0] = src_dimens[2] * src_strides[2];
00760 
00761     dst_strides[1] = element_size;
00762     dst_strides[2] = dst_dimens[1] * dst_strides[1];
00763     dst_strides[0] = dst_dimens[2] * dst_strides[2];
00764 
00765     weights = ComputeWeights(vpc, src_dimens[0], dst_dimens[0], filter_type);
00766     if (weights == NULL) {
00767         Dealloc(vpc, tmp1_array);
00768         Dealloc(vpc, tmp2_array);
00769         return(vpc->error_code);
00770     }
00771     code = vpResample(vpc, 3, src_dimens, dst_dimens, src_strides, dst_strides,
00772                       element_type, tmp2_array, out_array);
00773     Dealloc(vpc, weights);
00774 
00775     if (code != VP_OK) {
00776         Dealloc(vpc, tmp1_array);
00777         Dealloc(vpc, tmp2_array);
00778         return(code);
00779     }
00780 
00781     /* clean up */
00782     Dealloc(vpc, tmp1_array);
00783     Dealloc(vpc, tmp2_array);
00784     return(VP_OK);
00785 }
00786 
00787 /*
00788  * ComputeWeights
00789  *
00790  * Allocate and compute a filter weight table for a predefined filter type.
00791  */
00792 
00793 static float *
00794 ComputeWeights(vpc, src_xlen, dst_xlen, filter_type)
00795 vpContext *vpc;         /* context for storing table */
00796 int src_xlen;           /* number of samples in input scanline */
00797 int dst_xlen;           /* number of samples in output scanline */
00798 int filter_type;        /* type of filter (VP_BOX_FILTER, etc.) */
00799 {
00800     double scale_factor;
00801     int num_phases, num_taps, support, tap_limit, phases, table_size;
00802     int code;
00803     float *weights;
00804 
00805     switch (filter_type) {
00806     case VP_BOX_FILTER:
00807         support = 1;
00808         break;
00809     case VP_LINEAR_FILTER:
00810         support = 2;
00811         break;
00812     case VP_GAUSSIAN_FILTER:
00813         support = 3;
00814         break;
00815     case VP_BSPLINE_FILTER:
00816     case VP_MITCHELL_FILTER:
00817         support = 4;
00818         break;
00819     default:
00820         VPSetError(vpc, VPERROR_BAD_OPTION);
00821         return(NULL);
00822     }
00823     scale_factor = (double)dst_xlen / (double)src_xlen;
00824     if (scale_factor >= 1.0) {
00825         num_taps = support;
00826         num_phases = 1024;
00827     } else {
00828         num_taps = (double)support / scale_factor;
00829         tap_limit = 4;
00830         phases = 1024;
00831         while (1) {
00832             if (num_taps <= tap_limit) {
00833                 num_phases = phases;
00834                 break;
00835             }
00836             tap_limit *= 2;
00837             phases /= 2;
00838             if (phases <= 1) {
00839                 num_phases = 1;
00840                 break;
00841             }
00842         }
00843     }
00844     table_size = num_taps * num_phases * sizeof(float);
00845     Alloc(vpc, weights, float *, table_size, "weight_table");
00846     switch (filter_type) {
00847     case VP_BOX_FILTER:
00848         code = vpBoxFilter(num_taps, num_phases, weights, table_size);
00849         if (code != VP_OK) {
00850             Dealloc(vpc, weights);
00851             VPSetError(vpc, code);
00852             return(NULL);
00853         }
00854         break;
00855     case VP_LINEAR_FILTER:
00856         code = vpLinearFilter(num_taps, num_phases, weights, table_size);
00857         if (code != VP_OK) {
00858             Dealloc(vpc, weights);
00859             VPSetError(vpc, code);
00860             return(NULL);
00861         }
00862         break;
00863     case VP_GAUSSIAN_FILTER:
00864         code = vpGaussianFilter(VP_GAUSSIAN_SIGMA, num_taps, num_phases,
00865                                 weights, table_size);
00866         if (code != VP_OK) {
00867             Dealloc(vpc, weights);
00868             VPSetError(vpc, code);
00869             return(NULL);
00870         }
00871         break;
00872     case VP_BSPLINE_FILTER:
00873         code = vpBicubicFilter(1.0, 0.0, num_taps, num_phases, weights,
00874                                table_size);
00875         if (code != VP_OK) {
00876             Dealloc(vpc, weights);
00877             VPSetError(vpc, code);
00878             return(NULL);
00879         }
00880         break;
00881     case VP_MITCHELL_FILTER:
00882         code = vpBicubicFilter(1.0/3.0, 1.0/3.0, num_taps, num_phases,
00883                                weights, table_size);
00884         if (code != VP_OK) {
00885             Dealloc(vpc, weights);
00886             VPSetError(vpc, code);
00887             return(NULL);
00888         }
00889         break;
00890     }
00891     vpSetFilter(vpc, num_taps, num_phases, weights);
00892     return(weights);
00893 }
00894 
00895 /*
00896  * vpBoxFilter
00897  *
00898  * Compute filter weights for box filter.
00899  * For abs(x) < 0.5:
00900  *      k(x) = C
00901  * (C is chosen so that k(x) integrates to 1).
00902  * Otherwise:
00903  *      k(x) = 0
00904  */
00905 
00906 vpResult
00907 vpBoxFilter(num_taps, num_phases, weights, weights_bytes)
00908 int num_taps;   /* number of filter taps to compute */
00909 int num_phases; /* number of phases to compute */
00910 float *weights; /* array for storing filter weights
00911                    (num_taps*num_phases entries) */
00912 int weights_bytes; /* size of array (for error checking) */
00913 {
00914     int p, t;
00915     float *wptr;
00916     double value;
00917 
00918     if (weights_bytes != num_taps * num_phases * sizeof(float))
00919         return(VPERROR_BAD_SIZE);
00920     if (num_phases % 2 != 0)
00921         return(VPERROR_BAD_VALUE);
00922     wptr = weights;
00923     value = 1.0 / (double)num_taps;
00924     for (p = 0; p < num_phases; p++) {
00925         for (t = 0; t < num_taps; t++) {
00926             *wptr++ = value;
00927         }
00928     }
00929     return(VP_OK);
00930 }
00931 
00932 /*
00933  * vpLinearFilter
00934  *
00935  * Compute filter weights for linear interpolation.
00936  * For abs(x) < 1:
00937  *      k(x) = C * (1 - abs(x))
00938  * (C is chosen so that k(x) integrates to 1).
00939  * Otherwise:
00940  *      k(x) = 0
00941  */
00942 
00943 vpResult
00944 vpLinearFilter(num_taps, num_phases, weights, weights_bytes)
00945 int num_taps;   /* number of filter taps to compute */
00946 int num_phases; /* number of phases to compute */
00947 float *weights; /* array for storing filter weights
00948                    (num_taps*num_phases entries) */
00949 int weights_bytes; /* size of array (for error checking) */
00950 {
00951     int p, t;
00952     float *wptr1, *wptr2;
00953     double x0, delta_x, x, xa, tap_spacing, sum, normalize, value;
00954 
00955     if (weights_bytes != num_taps * num_phases * sizeof(float))
00956         return(VPERROR_BAD_SIZE);
00957     if (num_phases % 2 != 0)
00958         return(VPERROR_BAD_VALUE);
00959     wptr1 = weights;
00960     tap_spacing = 2.0 / (double)num_taps;
00961     x0 = -tap_spacing * ((double)num_taps/2.0 - 1.0);
00962     delta_x = tap_spacing / (double)num_phases;
00963     for (p = 0; p < num_phases/2; p++) {
00964         x = x0;
00965         sum = 0;
00966         for (t = 0; t < num_taps; t++) {
00967             if (x < 0.0)
00968                 xa = -x;
00969             else
00970                 xa = x;
00971             value = 1.0 - xa;
00972             wptr1[t] = value;
00973             sum += value;
00974             x += tap_spacing;
00975         }
00976         normalize = 1.0 / sum;
00977         for (t = 0; t < num_taps; t++) {
00978             wptr1[t] *= normalize;
00979         }
00980         wptr1 += num_taps;
00981         x0 -= delta_x;
00982     }
00983     wptr2 = wptr1;
00984     for (p = 0; p < num_phases/2; p++) {
00985         for (t = 0; t < num_taps; t++) {
00986             *wptr1++ = *--wptr2;
00987         }
00988     }
00989     return(VP_OK);
00990 }
00991 
00992 /*
00993  * vpBicubicFilter
00994  *
00995  * Compute filter weights for a Mitchell bicubic.
00996  *
00997  * See Mitchell, D.P. and Netravali, A.N., "Reconstruction filters in
00998  * computer graphics," Proc. SIGGRAPH '88 (Computer Graphics V22 N4),
00999  * p. 221-8.
01000  */
01001 
01002 vpResult
01003 vpBicubicFilter(b_value, c_value, num_taps, num_phases, weights, weights_bytes)
01004 double b_value; /* b in the filter kernel equation */
01005 double c_value; /* c in the filter kernel equation */
01006 int num_taps;   /* number of filter taps to compute */
01007 int num_phases; /* number of phases to compute */
01008 float *weights; /* array for storing filter weights
01009                    (num_taps*num_phases entries) */
01010 int weights_bytes; /* size of array (for error checking) */
01011 {
01012     int p, t;
01013     float *wptr1, *wptr2;
01014     double x0, delta_x, x, xa, tap_spacing, sum, normalize, value;
01015 
01016     if (weights_bytes != num_taps * num_phases * sizeof(float))
01017         return(VPERROR_BAD_SIZE);
01018     if (num_phases % 2 != 0)
01019         return(VPERROR_BAD_VALUE);
01020     wptr1 = weights;
01021     tap_spacing = 4.0 / (double)num_taps;
01022     x0 = -tap_spacing * ((double)num_taps/2.0 - 1.0);
01023     delta_x = tap_spacing / (double)num_phases;
01024     for (p = 0; p < num_phases/2; p++) {
01025         x = x0;
01026         sum = 0;
01027         for (t = 0; t < num_taps; t++) {
01028             if (x < 0.0)
01029                 xa = -x;
01030             else
01031                 xa = x;
01032             if (xa < 1.0) {
01033                 value = (((12. - 9.*b_value - 6.*c_value)*xa - 18. +
01034                           12.*b_value + 6.*c_value)*xa*xa + 6. -
01035                          2.*b_value) * 1./6.;
01036             } else {
01037                 value = ((((-b_value - 6.*c_value)*xa + 6.*b_value +
01038                            30.*c_value)*xa - 12.*b_value - 48.*c_value)*xa +
01039                          8.*b_value + 24.*c_value)* 1./6.;
01040             }
01041             wptr1[t] = value;
01042             sum += value;
01043             x += tap_spacing;
01044         }
01045         normalize = 1.0 / sum;
01046         for (t = 0; t < num_taps; t++) {
01047             wptr1[t] *= normalize;
01048         }
01049         wptr1 += num_taps;
01050         x0 -= delta_x;
01051     }
01052     wptr2 = wptr1;
01053     for (p = 0; p < num_phases/2; p++) {
01054         for (t = 0; t < num_taps; t++) {
01055             *wptr1++ = *--wptr2;
01056         }
01057     }
01058     return(VP_OK);
01059 }
01060 
01061 /*
01062  * vpGaussianFilter
01063  *
01064  * Compute filter weights for a Gaussian.
01065  * For abs(x) <= 1.0:
01066  *      k(x) = C * exp(-x*x/(2*sigma*sigma))
01067  * (C is chosen so that k(x) integrates to 1).
01068  * Otherwise:
01069  *      k(x) = 0
01070  */
01071 
01072 vpResult
01073 vpGaussianFilter(sigma, num_taps, num_phases, weights, weights_bytes)
01074 double sigma;   /* standard deviation */
01075 int num_taps;   /* number of filter taps to compute */
01076 int num_phases; /* number of phases to compute */
01077 float *weights; /* array for storing filter weights
01078                    (num_taps*num_phases entries) */
01079 int weights_bytes; /* size of array (for error checking) */
01080 {
01081     int p, t;
01082     float *wptr1, *wptr2;
01083     double x0, delta_x, x, tap_spacing, sigma2_inv, sum, normalize, value;
01084 
01085     if (weights_bytes != num_taps * num_phases * sizeof(float))
01086         return(VPERROR_BAD_SIZE);
01087     if (num_phases % 2 != 0)
01088         return(VPERROR_BAD_VALUE);
01089     wptr1 = weights;
01090     sigma2_inv = -1.0 / (2.0 * sigma * sigma);
01091     tap_spacing = 2.0 / (double)num_taps;
01092     x0 = -tap_spacing * ((double)num_taps/2.0 - 1.0);
01093     delta_x = tap_spacing / (double)num_phases;
01094     for (p = 0; p < num_phases/2; p++) {
01095         x = x0;
01096         sum = 0;
01097         for (t = 0; t < num_taps; t++) {
01098             value = exp(x*x*sigma2_inv);
01099             wptr1[t] = value;
01100             sum += value;
01101             x += tap_spacing;
01102         }
01103         normalize = 1.0 / sum;
01104         for (t = 0; t < num_taps; t++) {
01105             wptr1[t] *= normalize;
01106         }
01107         wptr1 += num_taps;
01108         x0 -= delta_x;
01109     }
01110     wptr2 = wptr1;
01111     for (p = 0; p < num_phases/2; p++) {
01112         for (t = 0; t < num_taps; t++) {
01113             *wptr1++ = *--wptr2;
01114         }
01115     }
01116     return(VP_OK);
01117 }