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
00032
00033
00034
00035
00036
00037
00038
00039 #include <math.h>
00040 #include <stddef.h>
00041 #include <stdio.h>
00042 #include <stdlib.h>
00043
00044
00045
00046 #include "bspline_interp.h"
00047
00048
00049
00050
00051
00052
00053 static inline void fComputeWeights(long SplineDegree, float w, float *Weight);
00054 static inline void dComputeWeights(long SplineDegree, double w, double *Weight);
00055 #define FLOOR(x) ((long)(x))
00056
00057
00058
00059
00060
00061
00062
00063
00064 extern float fInterpolatedValue
00065 (
00066 float *Bcoeff,
00067 long Width,
00068 long Height,
00069 float x,
00070 float y,
00071 long SplineDegree
00072 )
00073 {
00074 float *p;
00075 float xWeight[6], yWeight[6];
00076 float interpolated, w ;
00077 long xIndex[6], yIndex[6];
00078 long i, j, k, xi, yi;
00079
00080 if (SplineDegree & 1L) {
00081 xi = (long) FLOOR(x);
00082 yi = (long) FLOOR(y);
00083 }
00084 else {
00085 xi = (long)FLOOR(x+0.5f);
00086 yi = (long)FLOOR(y+0.5f);
00087 }
00088
00089
00090 fComputeWeights(SplineDegree, x - (float)xi, xWeight);
00091 fComputeWeights(SplineDegree, y - (float)yi, yWeight);
00092
00093 if (xi < SplineDegree || xi > (Width - SplineDegree) ||
00094 yi < SplineDegree || yi > (Height - SplineDegree) )
00095 {
00096
00097
00098
00099 long Width2 = 2L * Width - 2L, Height2 = 2L * Height - 2L;
00100
00101
00102
00103 if (SplineDegree & 1L) {
00104 i = (long)FLOOR(x) - SplineDegree / 2L;
00105 j = (long)FLOOR(y) - SplineDegree / 2L;
00106 for (k = 0L; k <= SplineDegree; k++) {
00107 xIndex[k] = i++;
00108 yIndex[k] = j++;
00109 }
00110 }
00111 else {
00112 i = (long)FLOOR(x + 0.5f) - SplineDegree / 2L;
00113 j = (long)FLOOR(y + 0.5f) - SplineDegree / 2L;
00114 for (k = 0L; k <= SplineDegree; k++) {
00115 xIndex[k] = i++;
00116 yIndex[k] = j++;
00117 }
00118 }
00119
00120
00121 for (k = 0L; k <= SplineDegree; k++) {
00122 xIndex[k] = (Width == 1L) ? (0L) : ((xIndex[k] < 0L) ?
00123 (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
00124 : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
00125 if (Width <= xIndex[k]) {
00126 xIndex[k] = Width2 - xIndex[k];
00127 }
00128 yIndex[k] = (Height == 1L) ? (0L) : ((yIndex[k] < 0L) ?
00129 (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
00130 : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
00131 if (Height <= yIndex[k]) {
00132 yIndex[k] = Height2 - yIndex[k];
00133 }
00134 }
00135
00136 interpolated = 0.0f;
00137 for (j = 0; j <= SplineDegree; j++) {
00138 p = Bcoeff + (ptrdiff_t)(yIndex[j] * Width);
00139 w = 0.0f;
00140 for (i = 0; i <= SplineDegree; i++) {
00141 w += xWeight[i] * p[xIndex[i]];
00142 }
00143 interpolated += yWeight[j] * w;
00144 }
00145 return(interpolated);
00146 }
00147 else
00148 {
00149
00150
00151
00152 p = &Bcoeff[yi*Width + xi];
00153 switch(SplineDegree)
00154 {
00155 case 2L:
00156 return ((xWeight[0] * p[-Width-1] + xWeight[1] * p[-Width] +
00157 xWeight[2] * p[-Width+1])
00158 * yWeight[0] +
00159 (xWeight[0] * p[ -1] + xWeight[1] * p[ 0] +
00160 xWeight[2] * p[ 1])
00161 * yWeight[1] +
00162 (xWeight[0] * p[ Width-1] + xWeight[1] * p[ Width] +
00163 xWeight[2] * p[ Width+1])
00164 * yWeight[2]);
00165 break;
00166 case 3L:
00167 return ((xWeight[0]*p[ -Width-1] + xWeight[1]*p[ -Width ] +
00168 xWeight[2]*p[ -Width+1] + xWeight[3]*p[ -Width+2])
00169 * yWeight[0] +
00170 (xWeight[0]*p[ -1] + xWeight[1]*p[ 0] +
00171 xWeight[2]*p[ +1] + xWeight[3]*p[ +2])
00172 * yWeight[1] +
00173 (xWeight[0]*p[ Width-1] + xWeight[1]*p[ Width ] +
00174 xWeight[2]*p[ Width+1] + xWeight[3]*p[ Width+2])
00175 * yWeight[2] +
00176 (xWeight[0]*p[2*Width-1] + xWeight[1]*p[2*Width ] +
00177 xWeight[2]*p[2*Width+1] + xWeight[3]*p[2*Width+2])
00178 * yWeight[3]);
00179 break;
00180 case 4L:
00181 return ((xWeight[0]*p[-2*Width-2] + xWeight[1]*p[-2*Width-1] +
00182 xWeight[2]*p[-2*Width ] + xWeight[3]*p[-2*Width+1] +
00183 xWeight[4]*p[-2*Width+2])
00184 * yWeight[0] +
00185 (xWeight[0]*p[ -Width-2] + xWeight[1]*p[ -Width-1] +
00186 xWeight[2]*p[ -Width ] + xWeight[3]*p[ -Width+1] +
00187 xWeight[4]*p[ -Width+2])
00188 * yWeight[1] +
00189 (xWeight[0]*p[ -2] + xWeight[1]*p[ -1] +
00190 xWeight[2]*p[ 0] + xWeight[3]*p[ 1] +
00191 xWeight[4]*p[ 2])
00192 * yWeight[2] +
00193 (xWeight[0]*p[ Width-2] + xWeight[1]*p[ Width-1] +
00194 xWeight[2]*p[ Width ] + xWeight[3]*p[ Width+1] +
00195 xWeight[4]*p[ Width+2])
00196 * yWeight[3] +
00197 (xWeight[0]*p[ 2*Width-2] + xWeight[1]*p[ 2*Width-1] +
00198 xWeight[2]*p[ 2*Width ] + xWeight[3]*p[ 2*Width+1] +
00199 xWeight[4]*p[ 2*Width+2])
00200 * yWeight[4]);
00201 break;
00202 case 5L:
00203 return ((xWeight[0]*p[-2*Width-2] + xWeight[1]*p[-2*Width-1] +
00204 xWeight[2]*p[-2*Width ] + xWeight[3]*p[-2*Width+1] +
00205 xWeight[4]*p[-2*Width+2] + xWeight[5]*p[-2*Width+3])
00206 * yWeight[0] +
00207 (xWeight[0]*p[ -Width-2] + xWeight[1]*p[ -Width-1] +
00208 xWeight[2]*p[ -Width ] + xWeight[3]*p[ -Width+1] +
00209 xWeight[4]*p[ -Width+2] + xWeight[5]*p[ -Width+3])
00210 * yWeight[1] +
00211 (xWeight[0]*p[ -2] + xWeight[1]*p[ -1] +
00212 xWeight[2]*p[ 0] + xWeight[3]*p[ 1] +
00213 xWeight[4]*p[ 2] + xWeight[5]*p[ 3])
00214 * yWeight[2] +
00215 (xWeight[0]*p[ Width-2] + xWeight[1]*p[ Width-1] +
00216 xWeight[2]*p[ Width ] + xWeight[3]*p[ Width+1] +
00217 xWeight[4]*p[ Width+2] + xWeight[5]*p[ Width+3])
00218 * yWeight[3] +
00219 (xWeight[0]*p[ 2*Width-2] + xWeight[1]*p[ 2*Width-1] +
00220 xWeight[2]*p[ 2*Width ] + xWeight[3]*p[ 2*Width+1] +
00221 xWeight[4]*p[ 2*Width+2] + xWeight[5]*p[ 2*Width+3])
00222 * yWeight[4] +
00223 (xWeight[0]*p[ 3*Width-2] + xWeight[1]*p[ 3*Width-1] +
00224 xWeight[2]*p[ 3*Width ] + xWeight[3]*p[ 3*Width+1] +
00225 xWeight[4]*p[ 3*Width+2] + xWeight[5]*p[ 3*Width+3])
00226 * yWeight[5]);
00227 break;
00228 default:
00229 return 0.0f;
00230 }
00231 }
00232 }
00233
00234
00235
00236
00237 extern void fAffine
00238 (
00239 float *Bcoeff,
00240 float *Image,
00241 long Width,
00242 long Height,
00243 float a11,
00244 float a12,
00245 float a21,
00246 float a22,
00247 float xShift,
00248 float yShift,
00249 long Masking,
00250 long SplineDegree
00251 )
00252 {
00253 float *p;
00254 float xWeight[6], yWeight[6];
00255 float interpolated, w;
00256 long xIndex[6], yIndex[6];
00257 long Width2 = 2L * Width - 2L, Height2 = 2L * Height - 2L;
00258 long i, j, k, xi, yi, x1,y1;
00259 float x, y, x0,y0;
00260
00261 if (SplineDegree < 2L || SplineDegree > 5L)
00262 {
00263 printf("Invalid spline degree: %ld\n", SplineDegree);
00264 return;
00265 }
00266
00267 x0 = xShift;
00268 y0 = yShift;
00269 for (y1 = 0L; y1 < Height; y1++)
00270 {
00271
00272
00273 x = x0;
00274 y = y0;
00275 for (x1 = 0L; x1 < Width; x1++)
00276 {
00277
00278
00279 if (Masking && ((x <= -0.5) || (((double)Width - 0.5) <= x)
00280 || (y <= -0.5) || (((double)Height - 0.5) <= y)))
00281 Image[y1*Width + x1] = 0.0F;
00282 else
00283 {
00284 if (SplineDegree & 1L) {
00285 xi = (long) (x);
00286 yi = (long) (y);
00287 }
00288 else {
00289 xi = (long)(x+0.5f);
00290 yi = (long)(y+0.5f);
00291 }
00292
00293
00294 fComputeWeights(SplineDegree, x - (float)xi, xWeight);
00295 fComputeWeights(SplineDegree, y - (float)yi, yWeight);
00296
00297 if (xi < SplineDegree || xi > (Width - SplineDegree) ||
00298 yi < SplineDegree || yi > (Height - SplineDegree) ) {
00299
00300 if (SplineDegree & 1L) {
00301 i = (long)(x) - SplineDegree / 2L;
00302 j = (long)(y) - SplineDegree / 2L;
00303 for (k = 0L; k <= SplineDegree; k++) {
00304 xIndex[k] = i++;
00305 yIndex[k] = j++;
00306 }
00307 }
00308 else {
00309 i = (long)(x + 0.5f) - SplineDegree / 2L;
00310 j = (long)(y + 0.5f) - SplineDegree / 2L;
00311 for (k = 0L; k <= SplineDegree; k++) {
00312 xIndex[k] = i++;
00313 yIndex[k] = j++;
00314 }
00315 }
00316
00317
00318
00319 for (k = 0L; k <= SplineDegree; k++) {
00320 xIndex[k] = (Width == 1L) ? (0L) : ((xIndex[k] < 0L) ?
00321 (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
00322 : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
00323 if (Width <= xIndex[k]) {
00324 xIndex[k] = Width2 - xIndex[k];
00325 }
00326 yIndex[k] = (Height == 1L) ? (0L) : ((yIndex[k] < 0L) ?
00327 (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
00328 : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
00329 if (Height <= yIndex[k]) {
00330 yIndex[k] = Height2 - yIndex[k];
00331 }
00332 }
00333
00334 interpolated = 0.0f;
00335 for (j = 0; j <= SplineDegree; j++) {
00336 p = Bcoeff + (ptrdiff_t)(yIndex[j] * Width);
00337 w = 0.0f;
00338 for (i = 0; i <= SplineDegree; i++) {
00339 w += xWeight[i] * p[xIndex[i]];
00340 }
00341 interpolated += yWeight[j] * w;
00342 }
00343 Image[y1*Width + x1] = interpolated;
00344 }
00345 else {
00346
00347 p = &Bcoeff[yi*Width + xi];
00348 switch(SplineDegree)
00349 {
00350 case 2L:
00351 Image[y1*Width + x1] = ((xWeight[0] * p[-Width-1] + xWeight[1] * p[-Width] +
00352 xWeight[2] * p[-Width+1])
00353 * yWeight[0] +
00354 (xWeight[0] * p[ -1] + xWeight[1] * p[ 0] +
00355 xWeight[2] * p[ 1])
00356 * yWeight[1] +
00357 (xWeight[0] * p[ Width-1] + xWeight[1] * p[ Width] +
00358 xWeight[2] * p[ Width+1])
00359 * yWeight[2]);
00360 break;
00361 case 3L:
00362 Image[y1*Width + x1] = ((xWeight[0]*p[ -Width-1] + xWeight[1]*p[ -Width ] +
00363 xWeight[2]*p[ -Width+1] + xWeight[3]*p[ -Width+2])
00364 * yWeight[0] +
00365 (xWeight[0]*p[ -1] + xWeight[1]*p[ 0] +
00366 xWeight[2]*p[ +1] + xWeight[3]*p[ +2])
00367 * yWeight[1] +
00368 (xWeight[0]*p[ Width-1] + xWeight[1]*p[ Width ] +
00369 xWeight[2]*p[ Width+1] + xWeight[3]*p[ Width+2])
00370 * yWeight[2] +
00371 (xWeight[0]*p[2*Width-1] + xWeight[1]*p[2*Width ] +
00372 xWeight[2]*p[2*Width+1] + xWeight[3]*p[2*Width+2])
00373 * yWeight[3]);
00374 break;
00375 case 4L:
00376 Image[y1*Width + x1] = ((xWeight[0]*p[-2*Width-2] + xWeight[1]*p[-2*Width-1] +
00377 xWeight[2]*p[-2*Width ] + xWeight[3]*p[-2*Width+1] +
00378 xWeight[4]*p[-2*Width+2])
00379 * yWeight[0] +
00380 (xWeight[0]*p[ -Width-2] + xWeight[1]*p[ -Width-1] +
00381 xWeight[2]*p[ -Width ] + xWeight[3]*p[ -Width+1] +
00382 xWeight[4]*p[ -Width+2])
00383 * yWeight[1] +
00384 (xWeight[0]*p[ -2] + xWeight[1]*p[ -1] +
00385 xWeight[2]*p[ 0] + xWeight[3]*p[ 1] +
00386 xWeight[4]*p[ 2])
00387 * yWeight[2] +
00388 (xWeight[0]*p[ Width-2] + xWeight[1]*p[ Width-1] +
00389 xWeight[2]*p[ Width ] + xWeight[3]*p[ Width+1] +
00390 xWeight[4]*p[ Width+2])
00391 * yWeight[3] +
00392 (xWeight[0]*p[ 2*Width-2] + xWeight[1]*p[ 2*Width-1] +
00393 xWeight[2]*p[ 2*Width ] + xWeight[3]*p[ 2*Width+1] +
00394 xWeight[4]*p[ 2*Width+2])
00395 * yWeight[4]);
00396 break;
00397 case 5L:
00398 Image[y1*Width + x1] = ((xWeight[0]*p[-2*Width-2] + xWeight[1]*p[-2*Width-1] +
00399 xWeight[2]*p[-2*Width ] + xWeight[3]*p[-2*Width+1] +
00400 xWeight[4]*p[-2*Width+2] + xWeight[5]*p[-2*Width+3])
00401 * yWeight[0] +
00402 (xWeight[0]*p[ -Width-2] + xWeight[1]*p[ -Width-1] +
00403 xWeight[2]*p[ -Width ] + xWeight[3]*p[ -Width+1] +
00404 xWeight[4]*p[ -Width+2] + xWeight[5]*p[ -Width+3])
00405 * yWeight[1] +
00406 (xWeight[0]*p[ -2] + xWeight[1]*p[ -1] +
00407 xWeight[2]*p[ 0] + xWeight[3]*p[ 1] +
00408 xWeight[4]*p[ 2] + xWeight[5]*p[ 3])
00409 * yWeight[2] +
00410 (xWeight[0]*p[ Width-2] + xWeight[1]*p[ Width-1] +
00411 xWeight[2]*p[ Width ] + xWeight[3]*p[ Width+1] +
00412 xWeight[4]*p[ Width+2] + xWeight[5]*p[ Width+3])
00413 * yWeight[3] +
00414 (xWeight[0]*p[ 2*Width-2] + xWeight[1]*p[ 2*Width-1] +
00415 xWeight[2]*p[ 2*Width ] + xWeight[3]*p[ 2*Width+1] +
00416 xWeight[4]*p[ 2*Width+2] + xWeight[5]*p[ 2*Width+3])
00417 * yWeight[4] +
00418 (xWeight[0]*p[ 3*Width-2] + xWeight[1]*p[ 3*Width-1] +
00419 xWeight[2]*p[ 3*Width ] + xWeight[3]*p[ 3*Width+1] +
00420 xWeight[4]*p[ 3*Width+2] + xWeight[5]*p[ 3*Width+3])
00421 * yWeight[5]);
00422 break;
00423 }
00424 }
00425 }
00426 x += a11;
00427 y += a21;
00428 }
00429 x0 += a12;
00430 y0 += a22;
00431 }
00432 }
00433
00434
00435
00436 static inline void fComputeWeights(long SplineDegree, float w, float *Weight)
00437 {
00438 float t,t0,t1,w2,w4;
00439
00440 switch (SplineDegree) {
00441 case 2L:
00442 Weight[1] = 3.0f / 4.0f - w * w;
00443 Weight[2] = (1.0f / 2.0f) * (w - Weight[1] + 1.0f);
00444 Weight[0] = 1.0f - Weight[1] - Weight[2];
00445 break;
00446 case 3L:
00447 Weight[3] = (1.0f / 6.0f) * w * w * w;
00448 Weight[0] = (1.0f / 6.0f) + (1.0f / 2.0f) * w * (w - 1.0f) - Weight[3];
00449 Weight[2] = w + Weight[0] - 2.0f * Weight[3];
00450 Weight[1] = 1.0f - Weight[0] - Weight[2] - Weight[3];
00451 break;
00452 case 4L:
00453 w2 = w * w;
00454 t = (1.0f / 6.0f) * w2;
00455 Weight[0] = 1.0f / 2.0f - w;
00456 Weight[0] *= Weight[0];
00457 Weight[0] *= (1.0f / 24.0f) * Weight[0];
00458 t0 = w * (t - 11.0f / 24.0f);
00459 t1 = 19.0f / 96.0f + w2 * (1.0f / 4.0f - t);
00460 Weight[1] = t1 + t0;
00461 Weight[3] = t1 - t0;
00462 Weight[4] = Weight[0] + t0 + (1.0f / 2.0f) * w;
00463 Weight[2] = 1.0f - Weight[0] - Weight[1] - Weight[3] - Weight[4];
00464 break;
00465 case 5L:
00466 w2 = w * w;
00467 Weight[5] = (1.0f / 120.0f) * w * w2 * w2;
00468 w2 -= w;
00469 w4 = w2 * w2;
00470 w -= 1.0f / 2.0f;
00471 t = w2 * (w2 - 3.0f);
00472 Weight[0] = (1.0f / 24.0f) * (1.0f / 5.0f + w2 + w4) - Weight[5];
00473 t0 = (1.0f / 24.0f) * (w2 * (w2 - 5.0f) + 46.0f / 5.0f);
00474 t1 = (-1.0f / 12.0f) * w * (t + 4.0f);
00475 Weight[2] = t0 + t1;
00476 Weight[3] = t0 - t1;
00477 t0 = (1.0f / 16.0f) * (9.0f / 5.0f - t);
00478 t1 = (1.0f / 24.0f) * w * (w4 - w2 - 5.0f);
00479 Weight[1] = t0 + t1;
00480 Weight[4] = t0 - t1;
00481 }
00482 }
00483
00484
00485
00486
00487
00488
00489 extern double dInterpolatedValue
00490 (
00491 double *Bcoeff,
00492 long Width,
00493 long Height,
00494 double x,
00495 double y,
00496 long SplineDegree
00497 )
00498 {
00499 double *p;
00500 double xWeight[6], yWeight[6];
00501 double interpolated, w ;
00502 long xIndex[6], yIndex[6];
00503 long i, j, k, xi, yi;
00504
00505 if (SplineDegree & 1L) {
00506 xi = (long) FLOOR(x);
00507 yi = (long) FLOOR(y);
00508 }
00509 else {
00510 xi = (long)FLOOR(x+0.5);
00511 yi = (long)FLOOR(y+0.5);
00512 }
00513
00514
00515 dComputeWeights(SplineDegree, x - (double)xi, xWeight);
00516 dComputeWeights(SplineDegree, y - (double)yi, yWeight);
00517
00518 if (xi < SplineDegree || xi > (Width - SplineDegree) ||
00519 yi < SplineDegree || yi > (Height - SplineDegree) )
00520 {
00521 long Width2 = 2L * Width - 2L, Height2 = 2L * Height - 2L;
00522
00523
00524
00525
00526
00527
00528 if (SplineDegree & 1L) {
00529 i = (long)FLOOR(x) - SplineDegree / 2L;
00530 j = (long)FLOOR(y) - SplineDegree / 2L;
00531 for (k = 0L; k <= SplineDegree; k++) {
00532 xIndex[k] = i++;
00533 yIndex[k] = j++;
00534 }
00535 }
00536 else {
00537 i = (long)FLOOR(x + 0.5) - SplineDegree / 2L;
00538 j = (long)FLOOR(y + 0.5) - SplineDegree / 2L;
00539 for (k = 0L; k <= SplineDegree; k++) {
00540 xIndex[k] = i++;
00541 yIndex[k] = j++;
00542 }
00543 }
00544
00545
00546 for (k = 0L; k <= SplineDegree; k++) {
00547 xIndex[k] = (Width == 1L) ? (0L) : ((xIndex[k] < 0L) ?
00548 (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
00549 : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
00550 if (Width <= xIndex[k]) {
00551 xIndex[k] = Width2 - xIndex[k];
00552 }
00553 yIndex[k] = (Height == 1L) ? (0L) : ((yIndex[k] < 0L) ?
00554 (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
00555 : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
00556 if (Height <= yIndex[k]) {
00557 yIndex[k] = Height2 - yIndex[k];
00558 }
00559 }
00560
00561 interpolated = 0.0;
00562 for (j = 0; j <= SplineDegree; j++) {
00563 p = Bcoeff + (ptrdiff_t)(yIndex[j] * Width);
00564 w = 0.0f;
00565 for (i = 0; i <= SplineDegree; i++) {
00566 w += xWeight[i] * p[xIndex[i]];
00567 }
00568 interpolated += yWeight[j] * w;
00569 }
00570 return(interpolated);
00571 }
00572 else
00573 {
00574
00575
00576
00577 p = &Bcoeff[yi*Width + xi];
00578 switch(SplineDegree)
00579 {
00580 case 2L:
00581 return ((xWeight[0] * p[-Width-1] + xWeight[1] * p[-Width] +
00582 xWeight[2] * p[-Width+1])
00583 * yWeight[0] +
00584 (xWeight[0] * p[ -1] + xWeight[1] * p[ 0] +
00585 xWeight[2] * p[ 1])
00586 * yWeight[1] +
00587 (xWeight[0] * p[ Width-1] + xWeight[1] * p[ Width] +
00588 xWeight[2] * p[ Width+1])
00589 * yWeight[2]);
00590 break;
00591 case 3L:
00592 return ((xWeight[0]*p[ -Width-1] + xWeight[1]*p[ -Width ] +
00593 xWeight[2]*p[ -Width+1] + xWeight[3]*p[ -Width+2])
00594 * yWeight[0] +
00595 (xWeight[0]*p[ -1] + xWeight[1]*p[ 0] +
00596 xWeight[2]*p[ +1] + xWeight[3]*p[ +2])
00597 * yWeight[1] +
00598 (xWeight[0]*p[ Width-1] + xWeight[1]*p[ Width ] +
00599 xWeight[2]*p[ Width+1] + xWeight[3]*p[ Width+2])
00600 * yWeight[2] +
00601 (xWeight[0]*p[2*Width-1] + xWeight[1]*p[2*Width ] +
00602 xWeight[2]*p[2*Width+1] + xWeight[3]*p[2*Width+2])
00603 * yWeight[3]);
00604 break;
00605 case 4L:
00606 return ((xWeight[0]*p[-2*Width-2] + xWeight[1]*p[-2*Width-1] +
00607 xWeight[2]*p[-2*Width ] + xWeight[3]*p[-2*Width+1] +
00608 xWeight[4]*p[-2*Width+2])
00609 * yWeight[0] +
00610 (xWeight[0]*p[ -Width-2] + xWeight[1]*p[ -Width-1] +
00611 xWeight[2]*p[ -Width ] + xWeight[3]*p[ -Width+1] +
00612 xWeight[4]*p[ -Width+2])
00613 * yWeight[1] +
00614 (xWeight[0]*p[ -2] + xWeight[1]*p[ -1] +
00615 xWeight[2]*p[ 0] + xWeight[3]*p[ 1] +
00616 xWeight[4]*p[ 2])
00617 * yWeight[2] +
00618 (xWeight[0]*p[ Width-2] + xWeight[1]*p[ Width-1] +
00619 xWeight[2]*p[ Width ] + xWeight[3]*p[ Width+1] +
00620 xWeight[4]*p[ Width+2])
00621 * yWeight[3] +
00622 (xWeight[0]*p[ 2*Width-2] + xWeight[1]*p[ 2*Width-1] +
00623 xWeight[2]*p[ 2*Width ] + xWeight[3]*p[ 2*Width+1] +
00624 xWeight[4]*p[ 2*Width+2])
00625 * yWeight[4]);
00626 break;
00627 case 5L:
00628 return ((xWeight[0]*p[-2*Width-2] + xWeight[1]*p[-2*Width-1] +
00629 xWeight[2]*p[-2*Width ] + xWeight[3]*p[-2*Width+1] +
00630 xWeight[4]*p[-2*Width+2] + xWeight[5]*p[-2*Width+3])
00631 * yWeight[0] +
00632 (xWeight[0]*p[ -Width-2] + xWeight[1]*p[ -Width-1] +
00633 xWeight[2]*p[ -Width ] + xWeight[3]*p[ -Width+1] +
00634 xWeight[4]*p[ -Width+2] + xWeight[5]*p[ -Width+3])
00635 * yWeight[1] +
00636 (xWeight[0]*p[ -2] + xWeight[1]*p[ -1] +
00637 xWeight[2]*p[ 0] + xWeight[3]*p[ 1] +
00638 xWeight[4]*p[ 2] + xWeight[5]*p[ 3])
00639 * yWeight[2] +
00640 (xWeight[0]*p[ Width-2] + xWeight[1]*p[ Width-1] +
00641 xWeight[2]*p[ Width ] + xWeight[3]*p[ Width+1] +
00642 xWeight[4]*p[ Width+2] + xWeight[5]*p[ Width+3])
00643 * yWeight[3] +
00644 (xWeight[0]*p[ 2*Width-2] + xWeight[1]*p[ 2*Width-1] +
00645 xWeight[2]*p[ 2*Width ] + xWeight[3]*p[ 2*Width+1] +
00646 xWeight[4]*p[ 2*Width+2] + xWeight[5]*p[ 2*Width+3])
00647 * yWeight[4] +
00648 (xWeight[0]*p[ 3*Width-2] + xWeight[1]*p[ 3*Width-1] +
00649 xWeight[2]*p[ 3*Width ] + xWeight[3]*p[ 3*Width+1] +
00650 xWeight[4]*p[ 3*Width+2] + xWeight[5]*p[ 3*Width+3])
00651 * yWeight[5]);
00652 break;
00653 default:
00654 return 0.0;
00655 }
00656 }
00657 }
00658
00659
00660
00661
00662 extern void dAffine
00663 (
00664 double *Bcoeff,
00665 double *Image,
00666 long Width,
00667 long Height,
00668 double a11,
00669 double a12,
00670 double a21,
00671 double a22,
00672 double xShift,
00673 double yShift,
00674 long Masking,
00675 long SplineDegree
00676 )
00677 {
00678 double *p;
00679 double xWeight[6], yWeight[6];
00680 double interpolated, w;
00681 long xIndex[6], yIndex[6];
00682 long Width2 = 2L * Width - 2L, Height2 = 2L * Height - 2L;
00683 long i, j, k, xi, yi, x1,y1;
00684 double x, y, x0,y0;
00685
00686 if (SplineDegree < 2L || SplineDegree > 5L)
00687 {
00688 printf("Invalid spline degree: %ld\n", SplineDegree);
00689 return;
00690 }
00691
00692 x0 = xShift;
00693 y0 = yShift;
00694 for (y1 = 0L; y1 < Height; y1++)
00695 {
00696 x = x0;
00697 y = y0;
00698 for (x1 = 0L; x1 < Width; x1++)
00699 {
00700 if (Masking && ((x <= -0.5) || (((double)Width - 0.5) <= x)
00701 || (y <= -0.5) || (((double)Height - 0.5) <= y)))
00702 Image[y1*Width + x1] = 0.0F;
00703 else
00704 {
00705
00706 if (SplineDegree & 1L) {
00707 xi = FLOOR(x);
00708 yi = FLOOR(y);
00709 }
00710 else {
00711 xi = FLOOR(x+0.5);
00712 yi = FLOOR(y+0.5);
00713 }
00714
00715
00716 dComputeWeights(SplineDegree, x - (double)xi, xWeight);
00717 dComputeWeights(SplineDegree, y - (double)yi, yWeight);
00718
00719 if (xi < SplineDegree || xi > (Width - SplineDegree) ||
00720 yi < SplineDegree || yi > (Height - SplineDegree) ) {
00721
00722 if (SplineDegree & 1L) {
00723 i = FLOOR(x) - SplineDegree / 2L;
00724 j = FLOOR(y) - SplineDegree / 2L;
00725 for (k = 0L; k <= SplineDegree; k++) {
00726 xIndex[k] = i++;
00727 yIndex[k] = j++;
00728 }
00729 }
00730 else {
00731 i = FLOOR(x + 0.5) - SplineDegree / 2L;
00732 j = FLOOR(y + 0.5) - SplineDegree / 2L;
00733 for (k = 0L; k <= SplineDegree; k++) {
00734 xIndex[k] = i++;
00735 yIndex[k] = j++;
00736 }
00737 }
00738
00739
00740
00741 for (k = 0L; k <= SplineDegree; k++) {
00742 xIndex[k] = (Width == 1L) ? (0L) : ((xIndex[k] < 0L) ?
00743 (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
00744 : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
00745 if (Width <= xIndex[k]) {
00746 xIndex[k] = Width2 - xIndex[k];
00747 }
00748 yIndex[k] = (Height == 1L) ? (0L) : ((yIndex[k] < 0L) ?
00749 (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
00750 : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
00751 if (Height <= yIndex[k]) {
00752 yIndex[k] = Height2 - yIndex[k];
00753 }
00754 }
00755
00756 interpolated = 0.0f;
00757 for (j = 0; j <= SplineDegree; j++) {
00758 p = Bcoeff + (ptrdiff_t)(yIndex[j] * Width);
00759 w = 0.0f;
00760 for (i = 0; i <= SplineDegree; i++) {
00761 w += xWeight[i] * p[xIndex[i]];
00762 }
00763 interpolated += yWeight[j] * w;
00764 }
00765 Image[y1*Width + x1] = interpolated;
00766 }
00767 else {
00768
00769 p = &Bcoeff[yi*Width + xi];
00770 switch(SplineDegree)
00771 {
00772 case 2L:
00773 Image[y1*Width + x1] =
00774 ((xWeight[0] * p[-Width-1] + xWeight[1] * p[-Width] +
00775 xWeight[2] * p[-Width+1])
00776 * yWeight[0] +
00777 (xWeight[0] * p[ -1] + xWeight[1] * p[ 0] +
00778 xWeight[2] * p[ 1])
00779 * yWeight[1] +
00780 (xWeight[0] * p[ Width-1] + xWeight[1] * p[ Width] +
00781 xWeight[2] * p[ Width+1])
00782 * yWeight[2]);
00783 break;
00784 case 3L:
00785 Image[y1*Width + x1] =
00786 ((xWeight[0]*p[ -Width-1] + xWeight[1]*p[ -Width ] +
00787 xWeight[2]*p[ -Width+1] + xWeight[3]*p[ -Width+2])
00788 * yWeight[0] +
00789 (xWeight[0]*p[ -1] + xWeight[1]*p[ 0] +
00790 xWeight[2]*p[ +1] + xWeight[3]*p[ +2])
00791 * yWeight[1] +
00792 (xWeight[0]*p[ Width-1] + xWeight[1]*p[ Width ] +
00793 xWeight[2]*p[ Width+1] + xWeight[3]*p[ Width+2])
00794 * yWeight[2] +
00795 (xWeight[0]*p[2*Width-1] + xWeight[1]*p[2*Width ] +
00796 xWeight[2]*p[2*Width+1] + xWeight[3]*p[2*Width+2])
00797 * yWeight[3]);
00798 break;
00799 case 4L:
00800 Image[y1*Width + x1] =
00801 ((xWeight[0]*p[-2*Width-2] + xWeight[1]*p[-2*Width-1] +
00802 xWeight[2]*p[-2*Width ] + xWeight[3]*p[-2*Width+1] +
00803 xWeight[4]*p[-2*Width+2])
00804 * yWeight[0] +
00805 (xWeight[0]*p[ -Width-2] + xWeight[1]*p[ -Width-1] +
00806 xWeight[2]*p[ -Width ] + xWeight[3]*p[ -Width+1] +
00807 xWeight[4]*p[ -Width+2])
00808 * yWeight[1] +
00809 (xWeight[0]*p[ -2] + xWeight[1]*p[ -1] +
00810 xWeight[2]*p[ 0] + xWeight[3]*p[ 1] +
00811 xWeight[4]*p[ 2])
00812 * yWeight[2] +
00813 (xWeight[0]*p[ Width-2] + xWeight[1]*p[ Width-1] +
00814 xWeight[2]*p[ Width ] + xWeight[3]*p[ Width+1] +
00815 xWeight[4]*p[ Width+2])
00816 * yWeight[3] +
00817 (xWeight[0]*p[ 2*Width-2] + xWeight[1]*p[ 2*Width-1] +
00818 xWeight[2]*p[ 2*Width ] + xWeight[3]*p[ 2*Width+1] +
00819 xWeight[4]*p[ 2*Width+2])
00820 * yWeight[4]);
00821 break;
00822 case 5L:
00823 Image[y1*Width + x1] =
00824 ((xWeight[0]*p[-2*Width-2] + xWeight[1]*p[-2*Width-1] +
00825 xWeight[2]*p[-2*Width ] + xWeight[3]*p[-2*Width+1] +
00826 xWeight[4]*p[-2*Width+2] + xWeight[5]*p[-2*Width+3])
00827 * yWeight[0] +
00828 (xWeight[0]*p[ -Width-2] + xWeight[1]*p[ -Width-1] +
00829 xWeight[2]*p[ -Width ] + xWeight[3]*p[ -Width+1] +
00830 xWeight[4]*p[ -Width+2] + xWeight[5]*p[ -Width+3])
00831 * yWeight[1] +
00832 (xWeight[0]*p[ -2] + xWeight[1]*p[ -1] +
00833 xWeight[2]*p[ 0] + xWeight[3]*p[ 1] +
00834 xWeight[4]*p[ 2] + xWeight[5]*p[ 3])
00835 * yWeight[2] +
00836 (xWeight[0]*p[ Width-2] + xWeight[1]*p[ Width-1] +
00837 xWeight[2]*p[ Width ] + xWeight[3]*p[ Width+1] +
00838 xWeight[4]*p[ Width+2] + xWeight[5]*p[ Width+3])
00839 * yWeight[3] +
00840 (xWeight[0]*p[ 2*Width-2] + xWeight[1]*p[ 2*Width-1] +
00841 xWeight[2]*p[ 2*Width ] + xWeight[3]*p[ 2*Width+1] +
00842 xWeight[4]*p[ 2*Width+2] + xWeight[5]*p[ 2*Width+3])
00843 * yWeight[4] +
00844 (xWeight[0]*p[ 3*Width-2] + xWeight[1]*p[ 3*Width-1] +
00845 xWeight[2]*p[ 3*Width ] + xWeight[3]*p[ 3*Width+1] +
00846 xWeight[4]*p[ 3*Width+2] + xWeight[5]*p[ 3*Width+3])
00847 * yWeight[5]);
00848 break;
00849 }
00850 }
00851 }
00852 x += a11;
00853 y += a21;
00854 }
00855 x0 += a12;
00856 y0 += a22;
00857 }
00858 }
00859
00860
00861 static inline void dComputeWeights(long SplineDegree, double w, double *Weight)
00862 {
00863 double t,t0,t1,w2,w4;
00864
00865 switch (SplineDegree) {
00866 case 2L:
00867 Weight[1] = 3.0 / 4.0 - w * w;
00868 Weight[2] = (1.0 / 2.0) * (w - Weight[1] + 1.0);
00869 Weight[0] = 1.0 - Weight[1] - Weight[2];
00870 break;
00871 case 3L:
00872 Weight[3] = (1.0 / 6.0) * w * w * w;
00873 Weight[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - Weight[3];
00874 Weight[2] = w + Weight[0] - 2.0 * Weight[3];
00875 Weight[1] = 1.0 - Weight[0] - Weight[2] - Weight[3];
00876 break;
00877 case 4L:
00878 w2 = w * w;
00879 t = (1.0 / 6.0) * w2;
00880 Weight[0] = 1.0 / 2.0 - w;
00881 Weight[0] *= Weight[0];
00882 Weight[0] *= (1.0 / 24.0) * Weight[0];
00883 t0 = w * (t - 11.0 / 24.0);
00884 t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
00885 Weight[1] = t1 + t0;
00886 Weight[3] = t1 - t0;
00887 Weight[4] = Weight[0] + t0 + (1.0 / 2.0) * w;
00888 Weight[2] = 1.0 - Weight[0] - Weight[1] - Weight[3] - Weight[4];
00889 break;
00890 case 5L:
00891 w2 = w * w;
00892 Weight[5] = (1.0 / 120.0) * w * w2 * w2;
00893 w2 -= w;
00894 w4 = w2 * w2;
00895 w -= 1.0 / 2.0;
00896 t = w2 * (w2 - 3.0);
00897 Weight[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - Weight[5];
00898 t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
00899 t1 = (-1.0 / 12.0) * w * (t + 4.0);
00900 Weight[2] = t0 + t1;
00901 Weight[3] = t0 - t1;
00902 t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
00903 t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
00904 Weight[1] = t0 + t1;
00905 Weight[4] = t0 - t1;
00906 }
00907 }
00908
00909
00910
00911
