# Digital image warping

Download 2.54 Mb.
 Page 28/30 Date 18.10.2016 Size 2.54 Mb.
 y hi 2(hi+h2) h2 y' hn-s 2(hn-s+hn_2) hn-2 6(r 1 - r0) 6(r 2 - r 1 ) 5(r._2 - r._3) The system of equations becomes determinable once the boundary conditions are specified. 289 A2.4.3. Boundary Conditions: Free-end, Cyclic, and Not-A-Knot A trivial choice for the boundary condition is achieved by setting y' =y'_ = O. This is known as thefree-end condition that results in natural spline interpolation. Since y' = 0, we know from Eq. (A2.2.6) that A2 = 0. As a result, we derive the following expression from Eq. (A2.3.1). Yl 3Ay0 y-t 2 - 2Ax (A2.4.5) Similarly, since Y-t = 0, 6A 3 + 2A2 = 0, and we derive the following expression from Eq. (A2.3.1). 2y_2 + 4y;-1 = 6 'Aye-2 (A2.4.6) AXn -2 Another condition is called the cyclic condition, where the derivatives at the end- points of the span are set equal to each other. y) = y_ (A2.4.7) y(; = y' The boundary condition that we shall consider is the not-a-knot condition. This requires y'" to be continuous across Xl and xn-2. In effect, this extrapolates the curve from the adjacent interior segments [de Boor 78]. As a result, we get A3 = A] (A2.4.8) x2 [ yo , ,] 1 +y] -2'x0 +Y+Yl j 'x2 [-Ayl ' 1 = _2_Xl +y 1 Replacing y with an expression in terms of y and y allows us to remain consistent with the stracture of a tridiagonal matrix already derived earlier. From Eq. (A2.4.2), we isolate y« and get lay0 Ayll ' AXl [ AXl+AX0] (A2.4.9) y : 3AXlL+J -yo c Substituting this expression into Eq. (A2.4.8) yields yAxt[Axo+AXt I +Y [AJ0+AXl] 2= AY [3AXoAXl +2Ax121 + AAYx [Ax 1 Similarly, the last row is derived to be 2 , A + -- [ 3AXn-3A Xn-2 + 2 Xn-31 A Xn-3 Xn-2 ß ß The version of this boundary condition expressed in terms of second derivatives is left to 290 the reader as an exercise. Thus far we have placed no restrictions on the spacing between the data points. Many simplifications are possible if we assume that the points are equispaced, i.e., Axk= 1. This is certainly the case for image reconstracfion, where cubic splines can be used to compute image values between regularly spaced samples. The not-a-knot boun- dary condition used in conjunction with the system of equations given in Eq. (A2.4.2) is shown below. To solve for the polynomial coefficients, the column vector containing the first derivatives must be solved and then substituted into Eq. (A2.3.1). 41 141 14 4 y y -Sy0 + 4Yl +Y2 3(y2 -Y0) 3(y3 -Yl) 3 (Vn - - y -3 ) --Yn-3 --4Yn-2 + 5yn-I (A2.4.10) A2.5. SOURCE CODE Below we include two C programs for interpolating cubic splines. The first pro- gram, called ispline, assumes that the supplied data points are equispaced. The second program, ispline_gen, addresses the more general case of irregularly spaced data. A2.5.1. Ispline The function ispline takes Y 1, a list of len 1 numbers in double-precision, and passes an interpolating cubic spline through that data. The spline is then resampled at len2 equal intervals and stored in list Y2. It begins by computing the ulknQwn first derivatives at each interval endpoint. It invokes the function getYD, which returns the first derivatives in the list YD. Along the way, function tridiag is called to solve the tridi- agohal system of equations shown in Eq. (A2.4.10). Since each derivative is coupled only to its adjacent neighbors on both sides, the equations can be solved in linear time, i.e., O(n). Once YD is initialized, it is used together with Y1 to compute the spline coefficients. In the interest of speed, the cubic polynomials are evaluated by using Horoer's rule for factoring polynomials. This requires three additions and three multipli- cations per evaluated point. Interpolating cubic spline function for equispaced points Input: Y1 is a list of equispaced data points with lenl entries Output: Y2 <- cubic spline sampled at len2 equispaced points ispline(Y1 ,lenl ,Y2,1en2) double*Y1, *Y2; int lenl, len2; int i, ip, oip; double *YD, A0, A1, A2, A3, x, p, fctr; /* compute 1st derivatives at each point -> YD */ YD = (double *) calloc(len1, sizeof(double)); getYD(Y1,YD,lenl); /* * p is real-valued position into spline * ip is interval's left endpoint (integer) * oip is left endpoint of last visited interval */ oip = -1; /* force coefficient initialization '/ fctr = (double) (lenl-1) / (len2-1); for(i=p=ip=0; i < len2; i++) { /* check if in new interval */ if(ip != oip) { /* update inte.'val */ oip = ip; /* compute spline coefficients */ A0 = Yl[ip]; A1 = YD{ip]; A2 = 3.0*(Yl[ip+l]-Yl[ip]) - 2.0*YD[ip] - YD[ip+l]; A3 = -2.0*(Yl[ip+l]-Yl[ip]) + YD[ip] + YD[ip+l]; } /* use Homer's rule to calculate cubic polynomial */ x=p-ip; Y2[i] = ((A3*x + A2)*x + A1)*x + A0; /* increment pointer */ ip = (p += fctr); } cfree((char *) YD); } 291 292 YD <- Computed 1st derivative of data in Y (len enidos) The not-a-knot boundary condition is used getYD(Y,YD,len) double *Y, *YD; int en; int i; YD = -5.0*Y + 4.0'Y + Y; for(i = 1; i < len-1; i++) Y D[i]=3.0*(Y[i+ 1 ]-Y[i- 1]); YD[len-1] = -Y[len-3] - 4.0*Y[len-2] + 5.0*Y[len-1]; /* solve for the tridiagonal matdx: YD=YD*inv(tridlag matrix) */ tridiag(YD,len); Linear time Gauss Elimination with backsubstitution for 141 tddiagonal matrix with column vector D. Result goes into D tridiag(D,len) double *D; int len; int i; double *C; ? init first two entries; C is dght band of tridiagonal */ C = (double *) calloc(len, sizeof(double)); D = 0.5*D; Dill = (Dill - D) / 2.0; C = 2.0; /* Gauss elimination; forward substitution */ for(i = 2; i < len-1; i++) { C[i] = 1.0 / (4.0 - C[i-1]); D[i] = (Dill - D[i-1]) / (4.0- C[i]); C[i] = 1.0/(4.0 - C[i-1]); D[i]= (Dill - 4*D[i-1]) / (2.0- 4*C[i]); /* backsubstitution */ for(i = len-2; i >= 0; i--) D[i] -= (D[i+l] * C[i+l]); cfree((char *) C); 293 A2.5.2. Ispline_gen The function ispline_gen takes the data points in (X1,Y1), two lists of len 1 numbers, and passes an interpolating cubic spline through that data. The spline is then resamplod at len 2 positions and stored in Y2. The resampling locations are given by X 2. The function assumes that X2 is monotonically increasing and lies withing the range of numbers in X 1. As beforo, we begin by computing the unknown first derivativos at each interval endpoint. The function getYD_gen is then invoked to return the first derivatives in the list YD. Along the way, function tridiagen is called to solve the tridiagonal system of equations given in Eq. (A2.4.2). Once YD is initialized, it is used together with Y 1 to compute the spline coefficients. Note that in this general case, additional consideration must now be givon to determine the polynomial interval in which the resampling point lios. Interpolating cubic spline function for irregularly-spaced points Input: Y1 is a list of irregular data points (lenl entries) Their x-coordinates are specified in Xl Output: Y2 <- cubic spline sampled according to X2 (len2 entries) Assume that Xl ,X2 entries are monotonically increasing isplinegen(X 1 ,Y1 ,lenl ,X2,Y2,1en2) double *Xl, *Y1, *X2, *Y2; inf lenl, len2; { int i, j; double *YD, A0, A1, A2, A3, x, dx, dy, pl, p2, p3; /* compute 1st derivatives at each point -> YD '/ YD = (double *) calloc(lenl, sizeof(double)); getYDgen(Xl ,Y1,YD,lenl); /* error checking '/ if(X2 < Xl II X2[len2-1] > Xl[lenl-1]) ( fprintf(stderr,"ispline_gen: Out of range0); exit(); } P * pl is left endpoint of interval * p2 is resampling position * p3 is dght endpoint of interval * j is input index for current interval */ p3 = X210] - 1; /* force coefficient initialization */ for(i=j=0; i < len2; i++) { /* check if in new interval "/ p2 = X2[i]; if(p2 > p3) { /* find the interval which contains p2 */ for(; jXl[j]; J++); if(p2 < Xl[I]) j--; pl = XI[J]; /* update left endpoint */ p3 = XI[j+I]; /* update right endpoint */ /* compute spline coefficients */ dx = 1.0/(Xl[j+I] - Xl[j]); dy = (Y1 [j+l] - YI[j]) * dx; A0: YI[j]; A1 = YD[j]; A2 = dx ' (3.0*dy - 2,0*YD[j] - YD[j+I]); A3 = dx*dx * (-2,0*dy + YD[j] + YD[j+I]); } /* use Homer's rule to calculate cubic polynomial '1 x=p2-pl; Y2[i] = ((A3*x + A2)*x + A1)*x + A0; } cfree((char *) ¾D); YD <- Computed 1st derivative of data in X,Y (len entdes) The not-a-knot boundary condition is used getYD_gen(X,Y,YD, len) double *X, *Y, *YD; int len; { int i; double h0, hl, r0, rl, *A, *B, *C; /* allocate memory for tridiagonal bands A,B,C */ A = (double *) calloc(len, sizeof(double)); B = (double *) calloc(len, slzeof(double)); C = (double ') calloc(len, sizeof(double)); /* init first row data */ h0 = X- X; hl = X- X; r0 = (Y - Y) / h0; rl = (Y - Y) / hl; B = hl * (h0+hl); C = (h0+hl) * (h0+hl); YD = r0*(3*h0*hl + 2*h1*h1) + rl*h0*h0; P init tddiagonal bands A, B, C, and column vector YD '/ /* YD will later be used to return the derivatives */ for(i = 1; i < len-1; i++) { h0 = X[i]- X[i-1]; hl = X[i+l] - X[i]; r0: (Y[i] - Y[i-1]) / h0; rl = (Y[i+l] - Y[i]) / hl; A[i] = hl; B[i]= 2 * (h0+hl); c[i] = h0; YD[i] = 3 * (r0*hl + rl*h0); I /* last row */ A[i] = (h0+hl) * (h0+hl); B[i]= h0 * (h0+hl); YD[i] = r0*hl*hl + rl*(3*h0*hl + 2*h0*h0); /* solve for the tridiagonal matrix: YD=YD*inv(tridiag matrix) */ tridiag_ge n(A,B,C,Y D,le n); cfree((char *) A); cfree((char *) B); dree((char *) C); Gauss Elimination with backsubstitution for general tridiagonal matrix with bands A,B,C and column vector D. 295 II 296 t ridiag_g en(A,B,C,D,len) double *A, *B, *C, *D; int len; int i; double b, *F; F = (double *) calloc(len, sizeof(double)); /* Gauss elimination; forward substitution */ b = B; D[O] = D[O] / b; for(i = 1;i < len; i++) { F[i] = C[i-1] / b; b = B[i] - Alii*Eli]; if(b == 0) { fpd ntf(stderr,"getY D_gen: divide-by-zero0); exit(); } D[i] = (D[i] - D[i-1]*A[i]) / b; ) f' backsubstitution */ for(i = len-2; i >= 0; i--) D[i] -= (D[i+l] * F[i+l]); cfme((char *) F); Appendix 3 FORWARD DIFFERENCE METHOD The method of forward differences is used to simplify the computation of polynomi- als. It basically extends the incremental evaluation, as used in scanline algorithms, to higher-order functions, e.g., quadratic and cubic polynomials. We find use for this in Chapter 7 where, for example, perspective mappings are approximated without costly division operations. In this appendix, we derive the forward differences for quadratic and cubic polynomials. The method is then generalized to polynomials of arbitrary degree. The (first) forward difference of a function f (x) is defined as Af(x) = f(x+)-f(x), >0 (A3.1) It is used together with f (x), the value at the current position, to determine f (x + ), the value at the next position. In our application, = 1, denoting a single pixel step size. We shall assume this value for in the discussion that follows. For simplicity, we begin by considering the forward difference method for linear interpolation. In this case, the first forward difference is simply the slope of the line passing through two supplied function values. That is, Af(x)=at for the function f(x)=a]x+ao. We have already seen it used in Section 7.2 for Gouraud shading, whereby the intensity value at position x+l along a scanline is computed by adding Af (x) to f (x). Surprisingly, this approach readily lends itself to higher-order interpo- lants. The only difference, however, is that Af (x) is itelf subject to update. That update is driven by a second increment, known as the second forward difference. The extent to which these increments are updated is based on the degree of the polynomial being evaluated. In general, a polynomial of degree N requires N forward differences. We now describe forward differencing for evaluating quadratic polyhomials of the form f(x) = a2x2 +alx +ao (A3.2) The first forward difference forf (x) is expressed as 297 298 Af(x) = f(x+l)--f(x) (A3.3) = a2(2x+ 1)+at Thus, Af (x) is a linear expression. If we apply forward differences to Af (x), we get A2f (x) = A(Af (x)) (A3.4) = Af(x+l)--Af(x) = 2a2 Since A2f (x) is a constant, there is no need for further terms. The second forward differ- ence is used at each iteration to update the first forward difference which, in turn, is added to the latest result to compute the new value. Each loop in the iteration can be rewritten as f(x+l) = f (x)+Af (x) (A3.5) A/(x+l) = Af (x)+ A2f (x) If computation begins at x = 0, then the basis for the iteration is given by f, A f, and A2f evaluated at x=0. Given these three values, the second-degree polynomial can be evaluated from 0 to lastx using the following C code. for(x = 0; x < lastx; x++) { f[x+l] = f[x] + Af; /* compute next point '/ Af += A2f; /* update 1st forward difference '/ } Notice that Afis subject to update by A2f, but the latter term remains constant throughout the iteration. A similar derivation is given for cubic polynomials. However, an additional for- ward difference constant must be incorporated due to the additional polynomial degree. For a third-degree polynomial of the form f(x) = a3x 3 +a2x 2 +alx +ao (A3.6) the first forward difference is 299 Af (x) = f (x + 1) --f (x) (A3.7) = 3a3x 2 +(3a3 +2a2)x +a3 +a2+al Since bf (x) is a second-degree polynomial, two more forward difference terms are derived. They are A2f(x) = z(Af(x)) = 6a3x+6a3 +2a2 (A3.8) A3f(x) = A(A2f(x)) = 6a 3 The use of forward differences to evaluate a cubic polynomial between 0 and lastx is demons'ated in the following C code. tor(x = 0; x < lastx; x++} { l[x+l] = f[2x] + zf; /* compute next point '/ Af += A f; /' update 1st forward difference*/ A2f += A3I; /* update 2nd forward difference */ } In contrast to the earlier example, this case has an additional forward difference term that must be updated, i.e., A2f. Even so, this method offers the benefit of computing a third-degree polynomial with only three additions per value. An alternate approach, using Homer's role for factoring polynomials, requires three additions and three multipli- cations. This makes forward differences the method of choice for evaluating polynomi- als. The forward difference approach for cubic polynomials is depicted in Fig. A3.1. The basis of the entire iteration is shown in the top row. For consistency with our discus- sion of this method in Chapter 7, texture coordinates are used as the function values. Thus, we begin with u0, AU0, and A2U0 defined for position x0, where the subscripts refer to the position along the scanline. In order to compute our next texture coordinate at x = 1, we add AU 0 to u 0. This is denoted by the arrows that are in contact with u0. Note that diagonal arrows denote addition, while vertical arrows refer to the computed sum, Therefore, u 1 is the result of adding AU 0 to u0. The following coordinate, u2, is the sum of ul and AU. The latter term is derived from zu0 and A2U0. This regular structure collapses nicely into a com- pact iteration, as demonstrated by the short programs given earlier. Higher-order polynomials are handled by adding more forward difference terms. This corresponds to augmenting Fig. A3.1 with additional columns to the right. The order of computation is from the left to right. That is, the summations corresponding to the diagonal arrows are executed beginning from the left coltann. This gives rise to the 300 adjacent elements directly below. Those elements are then combined in similar fashion. This cycle continues until the last diagonal is reached, denoting that the entire span of points has been evaluated. Position Value XO UO Xl Ul Au 0 X2 U 2 AU! A2U0 x3 u3 au2 a2u a3uo x4 u4 au3 2u2 3u Figure A3.1: Forward difference method. [Abdou 82] [Abram 85] [Akima 78] [Akima 84] [Anderson 90] [Andrews 76] [Antoniou 79] [Atteia 66] [Ballard 82] [Barnhill 77] REFERENCES Abdou, Ikram E. and Kwan Y. Wong, "Analysis of Linear Interpolation Schemes for Bi-Level Image Applications," IBM J. Res. Develop., vol. 26, no. 6, pp. 667-680, November 1982. Abram, Greg, Lee Westover, and Turner Whirred, "Efficient Alias-free Rendering Using Bit-Masks and Look-Up Tables," Computer Graphics, (SIGGRAPH '85 Proceedings), vol. i9, no. 3, pp. 53-59, July 1985. Aldma, H., "A Method of Bivariate Interpolation and Smooth Surface Fitting for Irregularly Distributed Data Points," ACM Trans. Math. Software, vol. 4, pp. 148-159, 1978. Akima, H., "On Estimating Partial Derivatives for Bivariate Interpola- tion of Scattered Data," Rocky Mountain J. Math., vol. 14, pp. 41-52, 1984. Anderson, Scott E. and Mark A.Z. Dippe, ' 'A Hybrid Approach to Facial Animation," ILM Technical Memo #1026, Computer Graphics Depart- ment, Lucasfilm Ltd., 1990. Andrews, Harry C. and Claude L. Patterson III, "Digital Interpolation of Discrete Images," IEEE Trans. Computers, vol. C-25, pp. 196-202, 1977. Antoniou, Andreas, Digital Filters: Analysis and Design, McGraw-Hill, New York, 1979. Atteia, M., "Existence et determination des foncfions splines a plusieurs variables," C. R. Acad. Sci. Paris, vol. 262, pp. 575-578, 1966. Ballard, Dana and Christopher Brown, Computer Vision, Prentice-Hall, Englewood Cliffs, NJ, 1982. Barnhill, Robert E., "Representation and Approximation of Surfaces," Mathematical Software III, Ed. by J.R. Rice, Academic Press, London, pp. 69-120, 1977. 301 302 REFmESCgS [Bennett 84a] [Bennett 84b] [Bergland 69] [Bernstein 71] [Bernstein 76] [Bier 86] [Bizais 83] [Blake 87] [Blinn 76] [BooIt 86a] [Boult 86b] [Braccini 80] [Bracewell 86] [Briggs 74] Bennett, Phillip P. and Steven A. Gabriel, "Spatial Transformation Sys- tem Including Key Signal Generator," U.S. Patent 4,463,372, Ampex Corp., July 3i, 1984. Bennett, Phillip P. and Steven A. Gabriel, "System for Spatially Transforming Images," U.S. Patent 4,472,732, Ampex Corp., Sep. 18, 1984. Bergland, Glenn D., "A Guided Tour of the Fast Fourier Transform," IEEE Spectrum, vol. 6, pp. 41-52, July i969. Bernstein, Ralph and Harry Silverman, "Digital Techniques for Earth Resource Image Data Processing," Proc. Amer. Inst. Aeronautics and Astronautics 8th Annu. Meeting, vol. C21, AIAA paper no. 71-978, October 1971. Bemstein, Ralph, "Digital Image Processing of Earth Observation Sen- sor Data," IBM J. Res. Develop., vol. 20, pp. 40-57, January 1976. Bier, Eric A. and Ken R. Sloan, Jr., "Two-Part Texture Mappings," IEEE Computer Graphics and Applications, vol. 6, no. 9, pp. 40-53, Sep- tember 1986. Bizais, Y., I.G. Zubal, R.W. Rowe, G.W. Bennett, and A.B. Brill, "2-D .Fitting and Interpolation Applied to Image Distortion Analysis," Pic- torial Data Analysis, Ed. by R.M. Haralick, NATO ASI Series, vol. F4, 1983, pp. 321-333. Blake, Andrew and Andrew Zisserman, Visual Reconstruction, MIT Press, Cambridge, MA, 1987. Blinn, James F. and Martin E. Newell, "Texture and Reflection in Com- puter Generated Images," Comm. ACM, vol. 19, no. 10, pp. 542-547, October 1976. Boult, Terrarice E., "Visual Surface Reconstraction Using Sparse Depth Data," Proc. IEEE Conference on Computer Vision and Pattern Recog- nition, pp. 68-76, June 1986. Boult, Terrance E., Information Based Complexity: Applications in Non- linear Equations and Computer Vision, Ph.D. Thesis, Dept. of Computer Sciense, Columbia University, NY, 1986. Braccini, Carlo and Giuseppe Marino, "Fast Geometrical Manipulations of Digital Images," Computer Graphics and Image Processing, vol. 13, pp. 127-141, 1980. Bracewell, Ron, The Fourier Transform and Its Applications, McGraw- Hill, NY, 1986. Briggs, I.C., "Machine Contouring Using Minimum Curvature," Geo- physics, vol. 39, pp. 39-48, 1974. [Brigham 88] [Burt 88] [Butler 87] [Caelli 81] [Casey 71] [Catmull 74] [Catmull 80] [Chen 88] [Clough 65] [Cochran 67] [Cook 84] [Cook 86] [Cooley 65] [Cooley 67a] [Cooley 67b] 303 Brigham, E. Oran, The Fast Fourier Transform and Its Applications, Prentice-Hall, Englewood Cliffs, NJ, 1988. Butt, Peter J., "Moment Images, Polynomial Fit Filters, and the Problem of Surface Interpolation," Proc. 1EEE Conference on Computer Vision and Pattern Recognition, pp. 144-152, June 1988. Butler, David A. and Patricia K. Pierson, "Correcting Distortion in Digi- Directory: filedownloadDownload 2.54 Mb.Share with your friends:

The database is protected by copyright ©ininet.org 2020
send message

Main page