00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "vp_global.h"
00032
00033
00034 #define FLTFRAC_TO_FIX31(f) ((int)((f) * 2147483648.))
00035
00036 typedef struct {
00037 int in_ptr_offset;
00038
00039
00040 float *wptr;
00041
00042 int tap_min;
00043 int tap_max;
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
00062
00063
00064
00065
00066 vpResult
00067 vpSetFilter(vpc, num_taps, num_phases, weights)
00068 vpContext *vpc;
00069 int num_taps;
00070 int num_phases;
00071 float *weights;
00072 {
00073 int num_ones, bit;
00074
00075
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
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
00095
00096
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;
00104 int *src_dimens;
00105 int *dst_dimens;
00106
00107 int *src_strides;
00108 int *dst_strides;
00109 int element_type;
00110
00111 void *in_array;
00112 void *out_array;
00113 {
00114 int num_taps;
00115 int num_phases;
00116 int in_x_count;
00117 int out_x_count;
00118 int in_x_stride;
00119 double scale_factor;
00120 double in_x0;
00121
00122 int index0;
00123
00124
00125 int phase0;
00126 int index_incr;
00127
00128 int phase_incr;
00129
00130 int unused_phase_bits;
00131
00132 FilterTemplate *template;
00133 float *weights;
00134 int in_offset;
00135 int index, phase;
00136 int out_x;
00137 int tap_min, tap_max;
00138 int bit, d;
00139
00140
00141 if (vpc->filter_weights == NULL)
00142 return(VPSetError(vpc, VPERROR_BAD_SIZE));
00143
00144
00145
00146
00147
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
00155
00156
00157
00158
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
00164
00165
00166
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
00187
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
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
00238
00239
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;
00247 int *src_dimens;
00248 int *dst_dimens;
00249
00250 int *src_strides;
00251 int *dst_strides;
00252 unsigned char *in_array;
00253 unsigned char *out_array;
00254 FilterTemplate *template;
00255 {
00256 int out_x;
00257 float *wptr;
00258 float acc;
00259 int tap;
00260 int tap_min, tap_max;
00261 unsigned char *in_ptr;
00262
00263 unsigned char *in_scan_ptr;
00264 unsigned char *out_ptr;
00265 unsigned char *out_scan_ptr;
00266 FilterTemplate *sample_template;
00267 int out_x_count;
00268 int in_x_stride;
00269 int out_x_stride;
00270 int *scan_coord;
00271 int done;
00272 int dim;
00273
00274
00275 out_x_count = dst_dimens[0];
00276 in_x_stride = src_strides[0];
00277 out_x_stride = dst_strides[0];
00278
00279
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
00286 in_scan_ptr = in_array;
00287 out_scan_ptr = out_array;
00288
00289 done = 0;
00290 while (!done) {
00291
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 }
00314
00315
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 }
00330
00331
00332 Dealloc(vpc, scan_coord);
00333 }
00334
00335
00336
00337
00338
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;
00346 int *src_dimens;
00347 int *dst_dimens;
00348
00349 int *src_strides;
00350 int *dst_strides;
00351 unsigned short *in_array;
00352 unsigned short *out_array;
00353 FilterTemplate *template;
00354 {
00355 int out_x;
00356 float *wptr;
00357 float acc;
00358 int tap;
00359 int tap_min, tap_max;
00360 unsigned short *in_ptr;
00361
00362 unsigned short *in_scan_ptr;
00363 unsigned short *out_ptr;
00364 unsigned short *out_scan_ptr;
00365 FilterTemplate *sample_template;
00366 int out_x_count;
00367 int in_x_stride;
00368 int out_x_stride;
00369 int *scan_coord;
00370 int done;
00371 int dim;
00372
00373
00374 out_x_count = dst_dimens[0];
00375 in_x_stride = src_strides[0];
00376 out_x_stride = dst_strides[0];
00377
00378
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
00385 in_scan_ptr = in_array;
00386 out_scan_ptr = out_array;
00387
00388 done = 0;
00389 while (!done) {
00390
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 }
00413
00414
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 }
00433
00434
00435 Dealloc(vpc, scan_coord);
00436 }
00437
00438
00439
00440
00441
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;
00449 int *src_dimens;
00450 int *dst_dimens;
00451
00452 int *src_strides;
00453 int *dst_strides;
00454 float *in_array;
00455 float *out_array;
00456 FilterTemplate *template;
00457 {
00458 int out_x;
00459 float *wptr;
00460 float acc;
00461 int tap;
00462 int tap_min, tap_max;
00463 float *in_ptr;
00464
00465 float *in_scan_ptr;
00466 float *out_ptr;
00467 float *out_scan_ptr;
00468 FilterTemplate *sample_template;
00469 int out_x_count;
00470 int in_x_stride;
00471 int out_x_stride;
00472 int *scan_coord;
00473 int done;
00474 int dim;
00475
00476
00477 out_x_count = dst_dimens[0];
00478 in_x_stride = src_strides[0];
00479 out_x_stride = dst_strides[0];
00480
00481
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
00488 in_scan_ptr = in_array;
00489 out_scan_ptr = out_array;
00490
00491 done = 0;
00492 while (!done) {
00493
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 }
00511
00512
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 }
00531
00532
00533 Dealloc(vpc, scan_coord);
00534 }
00535
00536
00537
00538
00539
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;
00546 int in_x, in_y;
00547 void *out_array;
00548 int out_x, out_y;
00549 int element_type;
00550
00551 int filter_type;
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
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
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
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
00633 Dealloc(vpc, tmp1_array);
00634 return(VP_OK);
00635 }
00636
00637
00638
00639
00640
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;
00647 int in_x, in_y, in_z;
00648 void *out_array;
00649 int out_x, out_y, out_z;
00650 int element_type;
00651
00652 int filter_type;
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
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
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
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
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
00782 Dealloc(vpc, tmp1_array);
00783 Dealloc(vpc, tmp2_array);
00784 return(VP_OK);
00785 }
00786
00787
00788
00789
00790
00791
00792
00793 static float *
00794 ComputeWeights(vpc, src_xlen, dst_xlen, filter_type)
00795 vpContext *vpc;
00796 int src_xlen;
00797 int dst_xlen;
00798 int filter_type;
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
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906 vpResult
00907 vpBoxFilter(num_taps, num_phases, weights, weights_bytes)
00908 int num_taps;
00909 int num_phases;
00910 float *weights;
00911
00912 int weights_bytes;
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
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943 vpResult
00944 vpLinearFilter(num_taps, num_phases, weights, weights_bytes)
00945 int num_taps;
00946 int num_phases;
00947 float *weights;
00948
00949 int weights_bytes;
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
00994
00995
00996
00997
00998
00999
01000
01001
01002 vpResult
01003 vpBicubicFilter(b_value, c_value, num_taps, num_phases, weights, weights_bytes)
01004 double b_value;
01005 double c_value;
01006 int num_taps;
01007 int num_phases;
01008 float *weights;
01009
01010 int weights_bytes;
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
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072 vpResult
01073 vpGaussianFilter(sigma, num_taps, num_phases, weights, weights_bytes)
01074 double sigma;
01075 int num_taps;
01076 int num_phases;
01077 float *weights;
01078
01079 int weights_bytes;
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 }