Interactive Multigrid Image Warping

German Cheung

Please contact German Cheung at for further details

Return to Homepage RealTime 3D Reconstruction Temporal SFS Human Kinematic Modeling and Motion Capture Human Motion Transfer


The major objective of the project is to utilize the fast nature of multigrid to image warping and build an interactive image warping program with both good local controls and smooth transition.


There are two major ways of performing image warping, namely mesh warping approach [10] and field morphing method [1]. For the first method, a mesh with predefined resolution is placed over the image and the user moves a subset of the mesh points to some desired positions. A warp is then computed based on all the mesh points. The advantage of this method is its fast computation speed. However, as the warping is calculated on a predefined resolution and each mesh point exert the same amount of influence on the warp, the expressiveness and flexibility of the method are poor. For the field morphing method, instead of imposing a mesh on the image, several feature points or lines and their movements are identified by the user. The warping is then determined by creating a "force field" on the image pixels. Although this method allows the user to express and control the exact warping of the feature points and lines, it is relatively slow as each image point in the image is affected by all the warping force fields generated by the features.

The main aim of this project is to combine the advantages of the above two methods to build a real-time image warping system with good local control ability and smooth transition. The project idea is discussed in Section 3. The image warping system is described in Section 4 and some results are included in Section 5. Finally a brief conclusion is given in Section 6.


The basic theory of this project combines the idea of mesh and feature based warping. The user specified the transitions of a set of feature points in the image. Each feature point transition imposes a constraint on the warping mesh and the constraints propagate smoothly to the whole image. A smooth warping mesh (with resolution equal to the image resolution) is then calculated based on the constraints and the pixels are warped according to the warping mesh.

In this section, we will first explain how the problem of producing a smooth warping mesh based on several constraints is indeed equivalent to the problem of reconstructing a tex2html_wrap_inline1142 D surface based on scattered data points. Then a fast multigrid algorithm is introduced to solve the problem.

From Image Warping to Surface Reconstruction

Figure 1 shows an image with some feature points and their desired transitions. The original positions of the feature points are represented by tex2html_wrap_inline1144 while the warped position are denoted by tex2html_wrap_inline1146 Mathematically, we can express the relation between tex2html_wrap_inline1148 as


Figure 1: Image with feature points and transitions.

From the equations, we can see that the x and y directions are indeed decoupled and hence can be considered separately. If we define the warping (displacement) functions to be tex2html_wrap_inline1150 and tex2html_wrap_inline1152 for the x and y directions respectively, the problem can be restated as finding smooth surfaces tex2html_wrap_inline1154 and tex2html_wrap_inline1156 such that


where N is the number of feature points. Although there exists many different methods for surface fitting [4], a fast method suggested by Terzopoulos in [9] is used in order to achieve real-time performance.

Before formulating the problem into its mathematical form, we have to define the notion of smoothness. One way of quantifying the smoothness of a tex2html_wrap_inline1164 D surface tex2html_wrap_inline1166 is to use the thin plate model. The surface is modeled as a deformable thin plate and the degree of smoothness is measured by the amount of strain energy tex2html_wrap_inline1168 . The lower the strain energy, the smoother is the surface. tex2html_wrap_inline1170 can be expressed as


where tex2html_wrap_inline1178 is the region of interest, tex2html_wrap_inline1180 and tex2html_wrap_inline1182 are the mean and Gaussian curvatures of the surface respectively while tex2html_wrap_inline1184 is the second derivative of tex2html_wrap_inline1186 with respect x.

By applying the model to our problem and expressing the constraints (Eqs (1)) as energy of fitting, the solution of the surface mesh tex2html_wrap_inline1188 (similarly for tex2html_wrap_inline1190 ) can be obtained by solving the optimization problem as follows [8] :


where the last term indicates how well the mesh is fitting the constraints with tex2html_wrap_inline1204 as the "stiffness" constant. Theoretically, the optimization problem can be solved by using variational calculus. However in that case, tex2html_wrap_inline1206 will be represented by the summation of some continuous basis functions [8] which is inappropriate for our application. Instead, since we are dealing with digital images, the discrete version of the optimization can be formulated by approximating the second derivative by finite differences between adjacent pixels. If we represent the discrete version of tex2html_wrap_inline1208 by tex2html_wrap_inline1210 , where tex2html_wrap_inline1212 and h is the grid size, the minimization becomes


where tex2html_wrap_inline1218 is the domain of interest and tex2html_wrap_inline1220 are the positions of the constraints in the discrete domain. Following the arguments in [8], the optimization problem is solved by finding the solution of the following linear equation







where tex2html_wrap_inline1222 equals 1 if there is a constraint at tex2html_wrap_inline1224 and 0 otherwise. The expressions for tex2html_wrap_inline1226 are computed in [8] and can be represented as the templates in Figure 2 below. Note that the templates are different for elements close to the boundary.

Figure 2: Templates for the matrix tex2html_wrap_inline1228 .


Although standard relaxation methods such as Gauss-Seidel can be used to solve the above linear equation, the convergence rate is too slow and is thus inappropriate in our interative application. Instead, the fast multigrid method can be used. As discussed in [3], for a linear problem, the multigrid method speeds up convergence by relaxing the residual equation at a lower resolution (coarser grid) and correcting the solution at the fine level as follows :

Correction Scheme

However, in our application, the above scheme will not work as the constraints are defined on the function tex2html_wrap_inline1258 but not on the error tex2html_wrap_inline1260 . This means that the the accuracy of the solution at a finer level cannot be consistently propagated to the coarser level by this correction scheme. As suggested in [7], this problem can be overcome by using a different multigrid scheme called full approximation scheme as first published in [2].

Full Approximation Scheme

Full Multigrid Full Approximation Scheme

The algorithm for the surface fitting can be summarized as follow :

  1. Propagate the constraints at the finest level tex2html_wrap_inline1296 down to the coarser levels tex2html_wrap_inline1298 .

  2. Solve the equations at the coarsest level tex2html_wrap_inline1300 directly or by iterative relaxation.

  3. At a level tex2html_wrap_inline1302 with k > 1 , obtain an initial guess from the coarser level tex2html_wrap_inline1306 by interpolation.

  4. Invoke the Full approximation Scheme at tex2html_wrap_inline1308 as discussed in 3.2.2.

  5. If k = L, stop, otherwise, set tex2html_wrap_inline1312 and goto Step 3.

The whole process is illustrated in Figure 3.

Figure 3: Illustration of the Full Multigrid Full Approximation Scheme.

The Image Warping System


An image warping system is built with an user-interface as shown in Figure 4. The left image is the original image for warping and the right image is the warped image. The user can arbitrarily click feature points and specify their warped positions. The transitions will be shown by colored arrows with square dots indicating the final positions of the feature points. The warping system can be operated in two modes, namely the fixed mode and the interactive mode. The fixed mode performs warping after the user has specified all the feature points and their transitions. The interactive mode operates when the user wants to see the warping interactively by dragging a specific feature point. A hard copy feature is also added to the system so that the user can output not only the warped image but also the original image with the transition arrows in TIFF format. There is also a clear button so that the user can clear all the marked feature points. An accuracy slider is provided as an interface for the user to control how accurate the warping functions should be. An example is also included in Figure 4.

Parameters setting

To prevent unpredictable warping, besides the constraints provided by the user, the four corners of the image are fixed. Moreover as in this application, it is desirable that the warping at the feature points are as close as possible to that specified by the user. Hence a large value of tex2html_wrap_inline1316 is chosen (note that tex2html_wrap_inline1318 is dependent on the grid size as discussed in [7]).

Figure 4: The image warping system.

Experimental Results

Surface Fitting

This part of the results illustrate the effectiveness of the multigrid method as compared to Gauss-Seidel. Two experiments are performed to compare the effectiveness of both methods for a (finest) gird size of 32 x 32. A constraint is imposed on the surface at the point (19, 16) with an amplitude of 10. This is shown in Figure 5(a). Figure 5(b) and Figure 5(c) show the results of using Gauss-Seidel method and Full Multigrid method respectively. Both methods have the same absolute errors of 138.428 Gauss Seidel takes 2.32s for the computation while Full Multigrid uses 0.27s to reach the same absolute error. We can compare these two surfaces with that in Figure 5(d) which is the correct solution for the problem. It is obvious that although the surfaces by Gauss-Seidel and Full Multigrid have the same absolute errors, the latter one is in fact much closer to the correct solution as compared to the former one. An additional constraint is added to the system and the tests are performed again. The results are shown in Figure 6 which indicates that the Full Multigrid grid method is superior over Gauss-Seidel method both in accuracy and speed.

Figure 5: 3-D plots of (a). the constraints, (b). approximate solutions by Gauss-Seidel, (c). Multigrid and (d). the correct solution.

Figure 6: 3-D plots of (a). the constraints, (b). approximate solutions by Gauss-Seidel, (c). Multigrid and (d). the correct solution.

Image Warping

For example, in order to find the warping function for one dimension (say x) for a 512 x 512 image by using Full Multigrid Full Approximation Scheme (with the fastest parameter setting), it takes about 3s on a SGI O2 with a 175 MHz IP32 processor. The time is reduced to 0.7s for a 256 x 256 image. Hence the update rate is about 1 Hz for the interactive mode. Real time warping is possible for 256x256 image if a fast enough machine is used.

A set of original and warped images are shown in Figure 7 as an indication of the quality of the warping. It can be seen that the warping is quit smooth among all the constraints.

 corig  cwarped
Figure 7: An example of warping a cat image using our warping system.


An interactive image warping system is built based on a surface fitting paradigm and fast multigrid method with smooth transitions and flexible local control. The system is able to perform close to real-time warping operations on images with sizes up to 256x256 with smooth transition. Moreover the technique can be extended to create animations [6] and perform three dimensional volume morphing [5].


T. Beier, S. Neely, ``Feature-based Image Metamorphosis, '' SIGGRAPH'92 Proceedings, 26(2), pg. 35-42, July, 1992.

A. Brandt, ``Multi-level adaptive techniques (MLAT) for partial differential equations : ideas and software,'' Mathematical Software III, Rice. J.R. (ed.), Academic Press, New York, pg. 277-318, 1979.

W.L. Briggs, ''A Multigrid Tutorial,'' SIAM, 1987.

S. J. Gortler, M. F. Cohen, ``Hierarchical and Variational Geometric Modeling with Wavelets, '' 1995 Symposium on Interactive 3D Graphics pg. 35-42, 1995.

A. Lerios, C. Garfinkle, M. Levoy, ``Feature-based Volume Metamorphosis, '' SIGGRAPH'95 Proceedings, August, 1995.

P. Litwinowicz, L. Williams, ``Animating Images with Drawings, '' SIGGRAPH'94 Proceedings, pg. 409-412, July, 1994.

D. Terzopoulos, ``Multi-level Reconstruction of Visual Surfaces: Variational Principles and Finite Element Representations,'' MIT A.I. Memo No. 671, April, 1982.

D. Terzopoulos, ``Multilevel computational processes for visual surface reconstruction,'' Computer Vision, Graphics, and Image Processing, 24, pg.52-96 1983.

D. Terzopoulos, ``The Computation of Visible-Surface Representations,'' IEEE Trans. PAMI, Vol. 10, No. 4, pg. 417-438, 1988.

G. Wolberg, ``Digital Image Warping'', IEEE Computer Society Press, 1990.