AU682425B2

AU682425B2 - Method and apparatus for constructing intermediate images for a depth image from stero images

Info

Publication number: AU682425B2
Application number: AU11680/95A
Authority: AU
Inventors: Sergei V. Fogel
Original assignee: Eastman Kodak Co
Current assignee: Eastman Kodak Co
Priority date: 1995-02-10
Filing date: 1995-02-10
Publication date: 1997-10-02
Anticipated expiration: 2015-02-10
Also published as: AU1168095A

Description

AUSTRALIA

Patents Act COMPLETE SPECIFICATION

(ORIGINAL)

Class Int. Class Application Number: Lodged: Complete Specification Lodged: Accepted: Published: Priority Related Art: a 0 4 a Name of Applicant: Eastman Kodak Company Actual Inventor(s): Sergei V. Fogel Address for Service: PHILLIPS ORMONDE FITZPATRICK Patent and Trade Mark Attorneys 367 Collins Street Melbourne 3000 AUSTRALIA Invention Title: METHOD AND APPARATUS FOR CONSTRUCTING INTERMEDIATE IMAGES FOR A DEPTH IMAGE FROM STEREO IMAGES Our Ref 400694 POF Code: 4703/4703 The following statement is a full description of this invention, including the best method of performing it known to applicant(s): -1- 19- 66,645 METHOD AND APPARATUS FOR CONSTRUCTING INTERMEDIATE IMAGES FOR A DEPTH IMAGE FROM STEREO IMAGES REFERENCE TO MICROFICHE APPENDIX A microfiche appendix of source code in the language having 63 total frames and one total microfiche. The microfiche appendix included herewith provides source code suitable for implementing an embodiment of the present invention on a Sun Microsystem Sparc computer running the Unix operating system.

10 BACKGROUND OF THE INVENTION Field of the Invention The present invention relates to a method and an apparatus for constructing a depth image from a pair of stereo images by interpolating intermediate images and then interlacing them into a single image.

Description Of The Related Art Producing images that have the perception of depth has traditionally been accomplished by photographic methods. Integral and lenticular photography have a long history of theoretical consideration and demonstration, but have only been met with limited commercial success. Many of the elementary concepts supporting integral and lenticular photography have been known for many years (see Takanori Okoshi, Three-Dimensional Imaging Techniques, Academic Press, New York, 1976; and G. Lippman, "E'preuves re'versibles, Photographics integrales," Comptes Rendus, 146, 446-451, March 2, 1908). The term integral photography refers to the composition of the overall image as an integration of a

"I

2 66,645 large number of microscopically small photographic image components. Each photographic image component is viewed through a separate small lens usually formed as part of a mosaic of identical spherically-curved surfaces embossed or otherwise formed onto the front surface of a plastic sheet of appropriate thickness. This sheet is subsequently bonded or held in close contact with the emulsion layer containing the photographic image components. Lenticular photography could be considered a special case of integral photography where the small lenses are formed as sections of cylinders running the full extent of the print area in the vertical direction. A recent commercial attempt at a form of lenticular photography is the Nimslo camera which is now being manufactured by a Hong Kong camera works and sold as a 15 Nishika camera. A sense of depth is clearly visible, however, the images resulting have limited depth realism and appear to jump as the print is rocked or the viewer's vantage relative to the print is changed.

An optical method of making lenticular photographs is described by Okoshi. A photographic camera is affixed to a carriage on a slide rail which allows it to be translated in a horizontal direction normal to the direction of the desired scene. A series of pictures is .taken where the camera is translated between subsequent exposures in equal increments from a central vantage point to later vantage points either side of the central vantage point. The distance that the lateral vantage points are displaced from the central vantage point is dependent on the maximum angle which the lenticular material can project photographic image components contained behind any given lenticule before it begins to project photographic image components contained behind an adjacent lenticule. It is not necessary to include a picture from the central vantage point, in which case the number of images will be even. If a picture from the central vantage point is included, the number of images will be odd. The sum of the total number of views contained between and including the lateral vantage points will determine the number of photographic components which eventually will be contained behind each lenticule.

3 66,645 The negatives resulting from each of these views are then placed in an enlarger equipped with a lens of the same focal length the camera lens. Since the camera had been moved laterally between successive exposures as previously described, the positions of the images in the original scene will be seen to translate laterally across the film format. Consequently, the position of the enlarged images from the negatives also move laterally with respect to the center of the enlarger's easel as successive negatives are placed in the film gate. An assemblage consisting of a sheet of photographic material oriented with it's emulsion in contact with the flat back side of a clear O plastic sheet of appropriate thickness having lenticules embossed or otherwise formed into its other side, is placed 15 on the enlarger easel with the lenticular side facing the senlarger lens. The position of this sheet on the easel is adjusted until the field of the central image is centered on the center of this assemblage, and an exposure of the information being projected out of the enlarger lens is made through the lenticules onto the photographic emulsion.

Subsequently, negatives from the successive exposures are placed in the film gate and the position of this assemblage is readjusted on the easel to reposition each respective view to the center of the assemblage, and additional exposures of the information being projected from the enlarger lens are be made. When all the views contained between the lateral vantages have been similarly exposed on the emulsion through the lenticular plastic sheet, the film sheet can be separated from the lenticula: plastic sheet and developed. If the aperture diameter of the enlarger lens is set to equal the amount of lateral shift between alternate views, the space behind each lenticule will be found to be exactly filled with photographic image components. The final step in this process is to reassemble the photographic film and the plastic sheet again with intimate contact between the emulsion layer and the flat side of the lenticular plastics heet and so positioned laterally, so that the long strips of adjacent images resulting from exposures through the cylindrical lenticules are again

I

4 66,645 positioned in a similar manner under the lenticules for viewing.

Ideally, an integral or lenticular photograph displays an infinite number of different angular views from each lenslet or lenticule. This is impossible since each angular view must have a corresponding small finite area of exposed emulsion or other hard copy media. Consequently, as an upper limit, the number of views must not exceed the resolution limit of the hard copy media. In the Nimslo print, the number of views behind each lens is limited to four views, two of which were considered left perspective views and the remaining two as right perspective views.

This is well below the resolution limit of conventional photographic emulsions and allows for only two options of S 15 stereoscopic viewing perspective as the viewer's head moves laterally.

The optical multiple projection method was also utilized in many experiments by Eastman Kodak researchers and engineers in the 1960's and 1970's and produced a lenticular photo displayed on the front cover of the 1969 Annual Report to Stockholders. This print had a large number of cameras taking alternate views of the scene to provide a smooth transition of perspectives for each lens. It is possible that as many as 21 different angular views were present and the result is much more effective.

This method of image recording is called an "indirect" technique because the final print recording is indirectly derived from a series of two-dimensional image recordings.

A more modern method of creating a lenticular print or depth image is to photographcally capture two or more spatially separated images of the same scene and digitize the images, or capture digital images electronically. The digitized images are used to drive a printer to print the images electronically onto a print or transparency media to which a lenticular faceplate is attached for viewing. This method of creating a depth image is typified by U.S. Patent 5,113,213 and U.S.

w One problem with the above-described methods -I pC-ii IPI 5 66,645 is that quite often the images are captured at a location distant from the printing site and if more images of the scene are necessary to provide a realistic sense of depth, the process of capturing, etc. must be completely performed again. Since images are digital and are processed for printing by a computer it is possible to solve this problem by using the computer to create the extra needed images using interpolation. However, conventional interpolation is not suitable because of problems related to object occlusion. A processfpr creatinq the extra needed images is described in U.S. S_ However, this approach does not produce the best intermediate image possible because an optimal solution is not always found during the search for the features of the intermediate image 15 due to the multidimensional searching and the requirement for image smoothing when a multidimensional search is being .performed.

perf:m A second problem with the above systems is that plural images must be stored between which the interpolated intermediate images are created.

What is needed is a system of stereo image interpolation which will produce multiple intermediate images between a pair of captured images such that the intermediate image quality is better because more optimal solutions to the search are found requiring less images W.C smoothing.

SUMMARY OF THE INVENTION It is an" object of the present inventien toproduce a depth image that is more appealing to thei wer.

It is another object of the prese, invention to create several intermediate images fro single pair of stereo images.

It is also an objt of the present invention to constrain the image se process in a horizontal direction to improve/i erpolated image quality and reduce the need for fi ng in image gaps with image smoothing processes.

RA; It is an additional object of the present \n to al n r-irety f the stereo imagGw in the- I- cF- Sa- According to one aspect of the invention there is provided a method, comprising the steps of: capturing a pair of stereo images of a scene separated in a horizontal direction; and creating an intermediate image between the pair by estimation of velocity vector fields constrained in the horizontal direction.

According to another aspect of the invention there is provided an apparatus for producing a depth image, comprising: a capture device capturing a pair of stereo images; a computer connected to said capture device and creating an intermediate image by horizontally constrained correspondence interpolation and creating the depth image from the pair and the intermediate image; and a depth image output device connected to said computer and outputting the depth image.

These together with other features and advantages which will be subsequently apparent, reside in the details of construction and operation as more fully hereinafter described and claimed, reference being had to the accompanying drawings forming a part hereof, wherein like numerals refer to like parts throughout. It is to be understood that the particularity of the detailed 20 description does not limit or supersede the generality of the foregoing description.

C AVINWORD ANDRCWSPECISII 1680CL DOC 66,645 -vcicaL *aireczLon in singLe operat.i--, The present invention accomplishes the a ve objects by performing an interpolation operation be een a pair of stereo images to produce plural interme 'ate images which are then used to produce a depth image The interpolation operation involves estimati the velocity vector field of the intermediate image The estimation for a particular intermediate image inv ves constraining the search for the correspondences btween the two images in the horizontal direction allowin the system to arrive at a more optimal solution which as a result does not require excessive gap filling y smoothing and which at the same time vertically ali ns the entirety of the images.

Tese together with other objects and 15 advantages ich will be subsequently apparent, reside in the deta s of construction and operation as more fully herei fter described and claimed, reference being had to t accompanying drawings forming a part hereof, wherein *nmna-a-is refor to likeo parts throughout.

BRIEF DESCRIPTION OF THE DRAWINGS Figure 1 depicts the hardware components of the present invention; Figure 2 illustrates two dimensional searching during correspondence processing; S 25 Figure 3 depicts one dimensional searching during correspondence processing; Figure 4 is a block diagram of the processes of the present invention; and Figure 5 illustrates disparity vector field processing.

DESCRIPTION OF THE PREFERRED EMBODIMENTS The present invention is 'esigned to obtain images intermediate to a pair of stereo ages so that a depth image, such as a lenticular or barrier image, can be created. The present invention, as illustrated in Figure 1, is implemented using conventional image processing equipment. The equipment includes one or more cameras ZRAL, that captures stereo images on film. The film images are r digitized by a conventional digitizer 12. An electronic 7 66,645 camera can of course substitute for the film camera and digitizer. The digitizer images are supplied to a computer 14 which performs the process discussed herein to create intermediate images. The images are interlaced and provided to a display device 16. The display device 16 can be a CRT with a lenticular o. b4arrier strip faceplate or it can be a print or transparency printer which produces a print to which a lenticular or barrier strip faceplate is attached.

To obtain the most appealing perception of depth it is desirable to produce a depth image with the maximum number of angular views. It is also desirable to reduce the number of optically recorded images to a single O pair of stereo images and to reconstruct the images corresponding to the missing angular views (the intermediate S 15 images) from this stereo pair. To reconstruct the intermediate images, first the correspondences between the stereo pair are established, then these correspondences are used to estimate the intermediate images through an interpolation process. The correspondence between the stereo pair establishes the apparent movement in the two images caused by differences in the position of the camera •lens. By considering the change in position of a point as a vector a means of estimating images as if a camera was positioned intermediate between the two actual cameras can be established. Not all points will have corresponding points, primarily due to occlusions, and allowances must be S. made to reconstruct images with these occurences. This vector field is commonly referred to as the optical flow (see Robot Vision, B.K.P. Horn, The MIT Press, 1991). One of the well-known characterizations of optical flow is that depth is inversely proportional to the magnitude of the optical flow vector. In this instance the zero parallax plane occurs in object space where the optical flow vectors have magnitude zero, points unaffected by the change of camera position, such as the background.

In creating a pixel 20 (See Figure 2) of an intermediate image 22 the first step is to find correspondences between the pair of images 24 and 26. To find the corresponding portions of pair of images a vector 8 66,645 or line 28 passing through the pixel 20 being created in the intermediate image is used to designate different areas on the images 24 and 26 are to be compared. Essentially the point of incidence on the images 24 and 26 can be systematically moved about two dimensional areas 30 and 32 searching for the areas that have the highest correlation in spectral and other image characteristics. Once these points are found the correspondences can be used to interpolate the value of the intermediate pixel 20. The correspondences between two images in two dimensions does not have a unique solution and, as a result, a search process as illustrated in Figure 2 may not produce the best solution. In addition, 4 because the solution is not unique it is necessary to perform smoothing of the image to reduce or eliminate noise 15 artifacts created in the intermediate image 22.

The process of establishing correspondence when stereo images are involved can be enhanced and the solution improved by recognizing that the two images are of the same scene. The only difference being the horizontal 20 displacement of the views or points from which the images were captured. Instead of searching a two dimensional area as in Figure 2, a one dimensional search along the horizontal axis 38 of the images is performed as illustrated in Figure 3.

Given this understanding of the concept of the present invention, the problem being solved by the present S• invention can be stated in more detail.

Given a pair of stereo images (images of a scene taken from different viewpoints) estimates of intermediate images (images taken from some intermediate viewpoints) must be found and the depth image constructed by interlacing the intermediate images. To estimate the intermediate images, the correspondences between the stere6 image pair must be established, then these correspondences used to estimate the intermediate images through an interpolation process. The axes of the scene coordinate system are defined as the horizontal axis, the vertical axis, and the depth axis. The selection of the viewpoints is restricted in such a way that the image planes of the

D-'

9 66,645 pair of stereo images coincide with the plane of the scene spanning the horizontal and the vertical axes. The cameras used in obtaining the stereo image pair are assumed to be positioned parallel to the plane of the scene spaning the horizontal and vertical axes, and the chenres in the camera positions are assumed to be strictly i' ;e ;vizontal direction. It is assumed that the horln-. ,l And the vertical axes of each image of the stereo image pair coincide with the horizontal and the vertical axes of the scene. Under these assumptions the correspondence between the images of the stereo image pair can be characterized by the intermediate disparity vector field with the constant V vertical component relating to the alignment of the images and with the variable horizontal component relating to the 15 depths of the corresponding points in the scene.

*To solve this problem as discussed hereinafter an overview of the right and the left images of the scene forming a stereo image pair and taken by a camera are digitized into multi-banded digital images with a fixed number of spectral bands, B. Each band of the right and the left digital images is defined as a two-dimensional array of numbers representing some measurements, known as pixels, taken from some grid of points n" on the sensed image. In the preferred embodiment, this grid is assumed to be a rectangular lattice. Given the right digital image 50 (See Figure 4) taken at the viewpoint tR and comprising of the spectral bands and given the left digital image 52 taken at the viewpoint tL and comprising of the spectral bands the process of constructing the depth image consisting of the bands is illustrated in Figure 4 and can be described as follows. First the correspondences between the left and right images are established 54, and the inte'mediate disparity vector fields V 1 (UK,VK) at some intermediate viewpoints tl,...,tK respectively are estimated 56 and 58. Then for each band b B and for each intermediate viewpoint the estimate of the intermediate image band Ib,k is obtained by the interpolation process. Finally, for each image band the depth image band Db is constructed by interlacing ~111 10 66,645 the intermediate image bands Ib,1 ,Ib,K to produce the depth image 62.

For each k 1, K the method of estimating the intermediate disparity vector field (UkVk) at the intermediate viewpoint tk is illustrated in Figure and can be described as the following multi-level resolution process: 1) Start 70 with the coarsest resolution level and the current estimate of the intermediate disparity vector field, (initially set to some constant); 2) Iteratively 72- 82 improve the current estimate; 3) Project 84 to the next finer level of resolution to obtain the initial estimate of the intermediate disparity vector field at the next finer level which is taken as the current estimate of the intermediate disparity vector field; 4) Repeat steps 2) and 15 3) until the finest resolution level is reached.

The proper estimate of the finest resolution level is taken as the final estimate of the intermediate disparity vector field.

On each level of the multi-level resolution '0 process a system of nonlinear equations relative to the unknown current estimate of the intermediate disparity vector field is formed 76, then iterative improvement of the current estimate of the intermediate disparity vector field is achieved by solving 76 this system of nonlinear equations. The system of nonlinear equations relative to the unknown current estimate of the intermediate disparity vector field is defined in terms of the unknown current estimate of the intermediate disparity vector field, the quantities related to the right and to the left pixels and the generalized spatial partial derivatives of the right and the left images. These quantities are obtained by filtering each band of the right and the left images with filters and with the spatial partial derivatives of these filters. The symbol G is used for the set of indices with each index g e G specifying a particular band and a particular filter. The multi-level resolution process (See U.S. application 631,750 incorporated by reference herein) is built around a multi-level resolution pyramid meaning that on each coarser resolution level the filtered righ and ~ara--- 11 66,645 the filtered left images, their spatial partial derivatives and the estimates of the intermediate disparity vector fields are defined on subgrids of grids at a finer resolution level.

At a given resolution level the bands of the filtered right images, the spatial partial derivatives of the bands of the filtered right images, the bands of the filtered left images, and the spatial partial derivatives of the bands of the filtered left images are defined on the same grid of points n' on the image plane of size N' by M'.

At the same resolution level the intermediate disparity vector fields are defined on a grid of points n on the image plane of size N by M. (In the following discussion the partial spatial derivatives taken with respect to the 15 vertical direction will have the subscript v, and analogously partial spatial derivatives with respect to the .:horizontal direction will have the subscript To form the system of nonlinear equations at this resolution level proceed as follows: For each index g G specifying a particular filter and a particular image band, form an optical flow function g, relative to a grid n as the difference between the estimated right filtered image band and the estimated left filtered image band. To obtain the estimated right filtered image band, find the grid of image points nR which are the taken from the viewpoint t R estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image p.ane from the viewpoint tk are the grid points n; and then interpolate the right filtered image fom the grid of points n' where the right filtered image is defined to the grid of points The estimated left filtered image band is determined analogously. To obtain the grid of image points each image point of the grid n is shifted by the amount equal to the scaled by the factor (tR-tk)/(tR-tL) value of the estimated intermediate disparity vector defined at that grid point. To obtain the grid of image points nL, each image point of the grid n is shifted by the amount equal to the scaled by the factor (tk-tL)/(tR-tL) value of the 12 66,645 estimated intermediate disparity vector defined at that grid point.

For each index g e G specifying a particular filter and a particular image band, form a spatial partial derivative g.u of the optical flow function gt with respect to the component u of the current estimate of the intermediate disparity vector corresponding to the index g as the defined on the grid n sum of the scaled by the factor (tR-tk)/(tR-tL) spatial rartial derivative of the estimated right filtered image band and the scaled by the factor (tktL)/(tR-tL) spatial partial derivative of the estimated left filtered image band. To obtain the spatial partial derivative of the estimated right filtered image band, find the grid of image points n, which are taken from the 15 viewpoint tR estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane from the viewpoint tk are the grid points n; and then interpolate the spatial partial derivatives of the right filtered image from the grid of points n' where the spatial partial derivative of the right filtered image is defined to the grid of points Sn,. To obtain the spatial partial derivative of the estimated left filtered image band, find the grid of image points n which are taken from the viewpoint tL estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane from the viewpoint t. are the grid points n; and then interpolate the spatial partial derivatives of the left filtered image from the grid of points n' where the spatial partial derivative of the left filtered image is defined to the grid of points n

L

To obtain the grid of image points nR, each image point of the grid n is shifted by the amount equal to the scaled by the factor (tR-tk)/(tR-tL) value of the estimated intermediate disparity vector defined at that grid point. To obtain the grid of image points n, each image point of the grid n is shifted by the amount equal to the scaled by the factor (tk-tL)/(tR-tL) value of the estimated intermediate disparity vector defined at that grid point.

II 13 66,645 For each index g E G specifying a particular filter and a particular image band, form a spatial partial derivative gt, of the optical flow function gt with respect to the component v of the current estimate of the intermediate disparity vector corresponding to the index g as the defined on the grid n sum of the scaled by the factor (tR-tk)/(tR-tL) spatial partial derivative with respect to the vertical component of the image coordinate system of the estimated right filtered image band corresponding to the index g and the scaled by the factor (tk-tL)/(tR-tL) spatial partial derivative with respect to the vertical component of the image coordinate system of the estimated left filtered image band corresponding to the index g. To obtain the spatial partial derivative ,ith respect to the vertical .s.e 15 component of the image coordinate system of the estimated oooo .'.".right filtered image band, find the grid of image points n

R

which are the taken from the viewpoint tR estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane from the viewpoint tk are the grid points a; and then interpolate the spatial partial derivatives with respect to the vertical component of the image coordinate system of the right filtered image from the grid of points S. n' where the spatial partial derivative with respect to the vertical component of the image coordinate system of the I right filtered image is defined to the grid of points RR- To obtain the spatial partial derivative with respect to the vertical component of the image coordinate system of the estimated left filtered image band, find the grid of image points OL which are taken from the viewpoint tL estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane from the viewpoint tk are the grid points n; and then interpolate the spatial partial derivatives with respect to the vertical component of the image coordinate system of the left filtered image from the grid of points n' where the spatial partial derivative with respect to the vertical component of the image coordinate system of the left filtered image is defined to the grid of points To L. 14 66,645 obtain the grid of image points nR, each image point of the grid n is shifted by the amount equal to the scaled by the factor (tR-tk)/(tR-tL) value of the estimated intermediate disparity vector defined at that grid point. To obtain the grid of image points each image point of the grid n is shifted by the amount equal to the scaled by the factor (tktL),/(tR-tL) value of the estimated intermediate disparity vector defined at that grid pioint.

Each non-boundary grid point of the rectangular grid n is surrounded by its eight nearest neighbors specifying eight different directions on the image plane. For each boundary grid point of the rectangular grid n the nearest neighbors are present only for some of these eight different directions. The symbol s is used to denote a vector on the image plane specifying one of these eight different directions on the image plane. The symbol S is used to denote the set of these eight vectors. For each vector s from the set S specifying a particular direction on the image plane, form a directional smoothness function 20 (s,Vu) for the variable horizontal component u of the S" current estimate of the intermediate disparity vector field as a finite difference approximation to the directional derivative of the horizontal component u of the current estimate of the intermediate disparity vector field. This directional smoothness function (s,Vu) is defined on the rectangilar subgrid f of the rectangular grid n with the property that every grid point of the rectangilar subgrid n8 has the nearest in the direction s neighbor on the rectangular grid n. At each grid point of the rectangular subgrid n. the directional smoothness function (s,Vu) is equal to the difference between the value of the horizontal component u of the current estimate of the intermediate disparity vector field at the nearest in the direction s neighbor of this grid point and the value of the horizontal component u of the current estimate of the intermediate disparity vector field at this grid point.

The system of nonlinear equations relative to the unknown current estimate of the intermediate disparity vector field is formed by combining the optical fjow

,II--I,

_111_ 15 66,645 function gt and its spatial partial derivatives gt gt for each filter and each image band specified by the index g e G together with the directional smoothness function Vu) for each image direction s e S and together with some constant parameters using four basic algebraic operations such as addition, subtraction, multiplication and division.

Given the right digital image taken at the right viewpoint tR and consistingof the bands and given the left digital image taken at the left viewpoint tL and consisting of the bands the method of estimating the intermediate digital image taken at the intermediate viewpoint tk and consisting of the bands II,k k based on the estimate of the intermediate disparity vector field (UkVk) obtained at the intermediate viewpoint tk can be described as follows. For each image band b B the estimate of the intermediate digital image band Ib, k is defined as the sum of the scaled by the factor (tk -tL)/(tR-tL) estimated right digital image band and the scaled by the factor (tR-tk)/(tR-tL) estimated left S 20 digital image band. To obtain the estimated right digital image band, find the grid of image points n R which are the taken from the viewpoint t R estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane from the viewpoint t k are the grid points f; and then interpolate the right digital image band from the grid of points n" where the right digital image is defined to the grid of points Rn. To obtain the estimated left digital image band, find the grid of image points nL which are the taken from the viewpoint tL estimated perspective projections onto the image plane of the visible points in the scene, whose perspective projections onto the image plane from the viewpoint tk are the grid points n; and then' interpolate the left digital image band from the grid of points n" where the left digital image is defined to the grid of points nL. To obtain the grid of image points nR, each image point of the grid n is shifted by the amount equal to the scaled by the factor (tR-tk)/(tR-tL) value of the estimated intermediate disparity vector defined at that I- rrrpl 16 66,645 grid point. To obtain the grid of image points nL, each image point of the grid n is shifted by the amount equal to the scaled by the factor (tk-tL)/(tR-tL) value of the estimated intermediate disparity vector defined at that grid point.

To provide a more thorough understanding of the fundamental operations of the present invention the following more detailed explanation is provided.

In what follows a more general problem is considered and some of the notation is changed. The above mentioned right digital image which is taken at the viewpoint tR and consists of the bands R, and the above mentioned left digital image which is taken at the Sviewpoint tL and consists of the bands are both represented in what follows by the irradiance image functions eS. For a fixed the above mentioned intermediate disparity vecotr field (Uk,Vk) which is defined at the intermediate viewpoint tk is represented in what follows by the parametric velocity vector field 20 estimate I En,ce[,oo]}. The pair of stereo viewpoints (tR, tL) is represented in what follows by the sequence of viewpoints where K"=l.

We shall deal with images of a three-dimensional scene and their spatial partial derivatives taken at the viewpoints within some interval T.

The initial images of a three-dimensional scene are, in general, discontinuous functions and their spatial partial derivatives are not defined at the points of discontinuity.

The points of discontinuity are, often, among the most important points in the images and cannot be easily ignored.

To overcome these difficulties the initial images are treated as generalized functions and their partial derivatives as generalized partial derivatives. These generalized functions are defined on some subset of the set of infinitely differentiable functions. The functions from the subset are called "testing functions". Parametric families of secondary images are then introduced through evaluations of the generalized functions associated with the initial images on the specific parametric family of 17 66,645 functions taken from the set of testing functions. The secondary images are infinitely differentiable, and their partial derivatives can be obtained through evaluations of the generalized functions associated with the initial images on the partial derivatives of this specific parametric family of functions. This process can be described as follows.

One way of constructing an initial image for a given viewpoint t e R (here R is a one-dimensional Euclidean space) is by projecting the light reflected by the objects in the scene onto a two-dimensional Euclidean space R 2 called the "image plane", and then identifying the irradiance value of each point in the image plane R 2 The function defined as above will be called the "irradiance image function". The value of the irradiance image function at the point e R 2 and the viewpoint t e R is assumed to be roughly proportional to the radiance of the point in the scene being imaged, which projects to such a point at the viewpoint t for every 20 e R 2 and every t e R. Different irradiance image

A

functions (x,y,t)e R 3 (here R 3 is a three-dimensional Euclidean space) can be obtained by changing the following aspects of the image formation process: the direction of a light source illuminating the scene, the color of such a light source, and the spectral responsivity function that is used to compute the irradiance.

Assume that each irradiance image function is locally integrable with respect to the Lebesgue measure dx dy dt in R 5 and thereby can be used to form a continuous linear functional rF (generalized function) defined on the locally convex linear topological space

#(R

3 The space $(R 3 consists of all infinitely differentiable functions having compact supports in the set

R

3 This means that for each function ~&e(R 3 there exists a closed bounded subset S eR 3 such that the function 0 is equal to zero at the points that are outside the subset S..

The topology of the space Q(R 3 is defined by a certain family of semi-norms. The functions p from the set $(R 3 1 18 66,645 will be called the "testing functions". The value of the generalized function F, associated with the irradiance image function t at the testing function Oet(R 3 is defined by the following relation: FP =ff f y, t)4) y, t) dx dy dt.

(2-1) The generalized functions r, associated with different irradiance image functions E(x,y,t) are united into the family (rEES). Here, to simplify the notation, the symbol O 10 E takes the role of the index (parameter) specifying the generalized function Fr associated with a particular irradiance image function in addition to its role as a part of the notation for such an irradiance image S* function. The symbol 3 denotes the set of indices.

S. 15 Another way of constructing an initi'al image at a given viewpoint t e R is by identifying th% .s.i of points M(t)cR 2 in the image plane, called the R'cstEure points", where significant variations in the projected light patterns represented by the irradiance image functions take 20 place at the viewpoint t, and then assigning feature value (x to each feature point The set of S. feature points M- UM(t) is assumed to be a closed subset teR Sof the space R 3 The function q(x,,y defined on the set M as above will be called the "feature image function".

Different feature image functions can be obtained by changing the criteria for selecting the set of feature points M and the criteria for assigning thefeature values The set M may, for example, be a finite combination of the following four types of the subsets of the space R 3 the three-dimensional regions, the two-dimensional regions, the ione-dimensional contours, and the isolated points.

By the family B B(M) of "Borel subsets" of

C-

19 66,645 the set M shall be meant the members of the smallest afinite, a-additive family of subsets of M which contains every compact set of M. Let be a a-finite, a-additive, and real-valued measure defined on the family B of Borel subsets of the set M. The feature image function (xyp,t) (xyptp) eM is assumed to be i-measurable and thereby can be used to form a continuous linear functional rF (generalized function) defined on the locally convex linear topological space D(R 3 The value of the generalized function F, associated with the feature image function r at the testing function pe$(R 3 is given by the following relation: (ow =fffr (xg ,y4, to) dH.

M

(2-2) 15 The generalized functions r associated with different feature image functions are united into the family {rjqeH}. Here, to simplify the notation, the symbol t takes the role of the index (parameter) specifying the 4 generalized function r, associated with a particular feature image function r7(x in addition to its role as a part of the notation for such a feature image function. The symbol H denotes the set of indices.

Given a generalized function F defined on the locally convex linear topological space #(R 3 and given 25 non-negative integer constants mx, my, mt the "generalized partial derivative" of the generalized axmxaymya t' function F is the generalized function F axmxayma am, defined on the locally convex linear topological space Q(R 3 as follows. The value of the generalized function I e I~i~ll~-ra~Psullasar 20 66,645 amx+my+m F at a testing function OE (R 3 is the value of axmxaymya tm, the generalized function F itself at the testing function ,m 3, (R 3 m y, t) +mx+mm 4Ym,+m y, t) y, t) ER 3 axxaaya"m A "combined generalized initial image function" will be defined as a linear combination of the generalized partial derivatives of the generalized functions associated with the irradiance image functions and of the generalized partial derivatives of the generalized functions associated with the feature image functions. Let A (A4, 1 n A,JES, i7eH) be a set of real-valued constants; let g(A) be an index attached to the set 1; and let mg(A){(m,x,my,m t, mrm,xm mr,tl ;|eE,yeH) be a set of non-negative integer constants corresponding to the set of constants 1. Then the combined generalized initial image 15 function corresponding to the set of constants A is the generalized function pg(A) defined on the locally convex linear topological space #(R 3 The value of the combined Sgeneralized initial image function corresponding to the set of constants I at the testing function Oe$(R 3 is given 20 by the relation I l l rt ws +E Ix c w (A axmx.ayyatn.( ax ay"y.atn.

(2-3) and is called the "observation" of the combined generalized initial image function corresponding to the set of constants A on the testing function 0. The combined generalized initial image functions corresponding to sets of constants A with different values are united into the family (rg(J)IAeA). Here the symbol A denotes the family of all different sets of constants A, while the symbol g r~P 21 66,645 denotes the one-to-one mapping from the family A onto the set of indices denoted by the symbol G. In what follows, to simplify the notation, the argument X a.pearing in the notation for a combined generalized initial image function is omitted and the symbol r 0 g E G is used instead of the symbol IeA to denote it.

Assume that there is a given fixed testing function *e#(R 3 which will be called the "measurement function". (An example of the measurement function is given later.) Let g e G; then for every point of some convex bounded subset n c R 2 (such points will be called "image points"), for every viewpoint t e T, and for every O value of the parameter ae[l,o), the value of the component of the image corresponding to the parameter value a at the image point and the viewpoint t is determined as the observation of the combined generalized initial image function rF on the testing function iOx,y,tE (R3) S 20 (2-4) Note that, to simplify the notation, the symbol g is used as Sthe part of the notation go for the component of the image corresponding to the parameter value a in addition to its role as the index for such a component.

25 The parametric vector valued function g9(x,y,t)=(go(x,y,t) geG}, (x,y,t)EnxT, e[l,oo) with each component go(x,y,t), geG defined on the se. nxT as the observation rg(~x,yt) of the combined generalized initial image function Fg on the testing function J1axy,tE$(R 3 specified by the relation for every parameter value a from the set (here the symbol n x T denotes the Cartesian product of the set n on the set T) will be called the "parametric viewpoint-varying image function". The points from the set n x T will be called "viewpoint-varying image points". By the "image" r sp -s sW~ llsRBBrs~-R- 22 66,645 corresponding to the parameter value ae[l,oo) and taken at the viewpoint t e T is meant the collection, (ga(x,y,t) of the values of the parametric viewpoint-varying image function ga(x,y,t), (x,y,t)enxT.

The collection (go(x,y,t) I (x,y)en,oa[l,w) of the images (x,y)en) corresponding to every parameter value e[l,co) and taken at a fixed viewpoint teT will be called the "parametric image" taken at the viewpo'nt t.

Each component gO geG of the parametric viewpoint-varying image function ga(x,y,t) is infinitely differentiable everywhere in the domain n x T, and its partial derivatives with respect to the variables x, y, t can be obtained as observations of the combined generalized initial image function rF on the partial derivatives with respect to the parameters x, y, t of the .r testing function *OytE$(R 3 specified by the relation For example, the components gy°(x,y,t), geG of the first-order partial derivatives IgeG}, gya(x,y,t)={gyo(x,y,t) jgeG) of the 20 parametric viewpoint-varying image function are given as the observations g (QW, g(ay c) of t combined generalized initial image function rF on the ao a.

testing functions ,y ,y.t (R3) for every En;<T, ae[l,o), where a a, 1 y /o, o) t)eR 3 a 1 y, V E) 1 tE)ER (2-6) II Ir I II m( 23 66, 645 and y,t) ilI y, t) are the partial derivatives of the function with respect to the variables x,y.

Let t-At-,t,t+At~cR be three consecutive viewpoint (here At' and At' are positive viewpoint increments), and let (x(t-At) y(t-At'))cR 2 be a projection of some point in the scene taken at the viewpoint t-At", while (x(t+At+),y(t+At+))6R 2 is a projection of the same point taken at the viewpoint t+At'. The vector (UA(XlyltAt,'At+),vA(x,y,t,At",At+)) defined by the relations uS(x",y (2-7) vA(x,,t,At,At)=(y(t+At) (2-8) where 15 X= (At x t +A t A x t-At (A t+ A t I v= (A t-Y t+A t At"Y t A t (At 29 :(2-10) will be called the "average velocity vector" of the scene projection at the point y) cR 2 and the viewpoint teR based on the scene projections at the viewpoint t+At+ and the viewpoint t-At-. The vector (uA y, t,At',At+) can be defined only at those points (x,y) from the set R 2 that are projections of the points in the scene visible at the viewpoint t-At' as well as the viewpoint t+At+. The collection of such image points will be denoted by the symbol fn(t,At',At 4 The "instantaneous velocity vector" of the scene projection at the point eR 2 and the viewpoint teR will now be defined. Let At-, At' be positive viewpoint increments, and let W(x,y,t,At*,At+) be the set which is equal to a single-element set containing the average velocity vector ~prrCB -r~u 24 66,645 (uA(X,y,t,At,At ),vA(x,y,t,AtAt if (x,y)en(tAt-,At and is equal to an empty set 0 otherwise. Then for every positive viewpoint increment At<6t (here St is a positive real constant), the set W(x,y,t,At) will be defined from the relation W(x, U W(x,y, t,A-,At), o<At-<At 0<At*<At (2-11) while the set W(x,y,t) will be defined as W(x,y, t) =nW(x,y, At) 0<At<t 10 (2-12) where the symbol W(x,y, t, At) means the topological closure of the set W(x,y,t,At). It will be assumed that W(x,y,t) are single-element sets for almost every (x,y,t)eR 3 This assumption means that the subset of the 15 set R 3 containing the points where the sets W(x,y,t) are not sinlte-element sets, has a Lebesgue measure equal to zero. For every point (x,y,t)eR 3 such that the set W(x,y,t) consists of a single element, that element is selected as the value of the instantaneous velocity vector 20 otherwise a zero vector is chosen for the value of the instantaneous velocity vector (ui(x,y,t) Assume that the functions ui(x,y,t) ,v 1 (x,y,t)eR 3 are locally integrable with respect to the Lebesgue measure dx dy dt in R 3 and thereby can be used to form continuous linear functionals U i

V

I

(generalized functions) defined on the locally convex linear topological space #(R 3 The values of the generalized functions U1, V I at the testing function e$ (R 3 are given by the relations I pO~q l~r~ II qllBd 25 5 -66f645 U, 4)=fffU, (X -V,0 Y(02-13)t

,R

3 (2-14) Let be an irnage point from 12, t a viewpoint from T, end a a parameter value from then the "velocity vector" of the image corresponding to the parazteter value a at the point and the viewpoint t is defined from the relations u(X, Y,L 0 lfffu(, 1 E4~-)a a3. R 3 crxdS' dE, (2-16) 5 3 fv,t is idpndn of the-x imag point were and is the marment fauncteio from the setle the 3 "vulociy t)v 0 fi inendt of the image poiyt n correspandin therefoe willabedtedalea n vtakwhile the collecin of The volelciyvctr (u7(m,y,t) vO(t)) I(xfy) eflaE(,o) of the velocity vector 26 66,645 fields I(x,y) n) of the images (x,y)en) corresponding to every parameter value e[l,co) and taken at a fixed viewpoint teT will be called the "parametric velocity vector field" of the parametric image (go(x,y,t) I en,oe taken at the viewpoint t.

Let to" t, t k be a finite increasing sequence of viewpoints from some viewpoint interval T, and let t be a viewpoint from the same viewpoint interval T.

Then the estimation problem can be formulated as follows.

Given a sequence of parametric images (go(X,y,tk/) I (x,y)Een,ae[l,) )k=K taken at the viewpoints tk, k=0, which will be called the "parametric viewpoint-varying image sequence", and given the viewpoint t, find the parametric vector field I En,E[l,oo) which is an estimate of the parametric velocity vector field I(x,y)en,aE[l,co)) of the parametric image (go(x,y,t) I (x,y)en,ae[l,co)) taken at the viewpoint t.

Although changes in the scene giving rise to S 20 the parametric velocity vector field (x,y)en,aet[l,co)} are reflected in the parametric viewpoint-varying image sequence (go(x,y,tk") I noral,o)k 0

K

the relation jetween them is not necessarily unique. The same parametric S 25 viewpoint-varying image sequence can be associated with different parametric velocity vector fields and, vice versa, the same parametric velocity vector field can result from different parametric viewpoint-varying image sequences.

The method of the present invention, then, includes determining an estimate of the parametric velocity vector field corresponding to a given parametric viewpoint-varying image sequence. In order for this determination to be possible, specific assumptions have to be made about the scene being imaged and about the imaging process itself. The assumptions we make are described in the following. Based on these assumptions constraints are imposed on the estimate of the parametric velocity vector field. The determination of the estimate is then reduced to solving the system of equations arising from such P Y A ~B L- 27 66,645 constraints for a given parametric viewpoint-varying image sequence.

Let (go(x,y,tk') en,aeCl,o))= 0 K" be a given parametric viewpoint-varying image sequence taken at the viewpoints tk', k=0, which form an increasing sequence within the viewpoint interval T, and let t be a given viewpoint within the viewpoint interval T. We consider the problem of finding an estimate I (x,y)en,ae[l,o)) of the parametric velocity vector field I en,aE(l,) as a function of the parametric viewpoint-varying image sequence (Xy) En,' [1l,o) In what follows, the estimate a I En,ae[l,o) of the parametric velocity vector field I (x,y)en,ae~[l,o) is determined as the parametric vector field satisfying a set of constraints. These constraints are based on the following assumptions: 1. The scene to be imaged is assumed to have near-constant and near-uniform illumination, which means that the changes in the incident illumination of each surface patch in the scene are small across the space, and are mainly due to the orientation of such surface patch relative to the light source.

2. The radiance of each surface point in the scene is assumed to be nearly Lambertian, which means that locally its dependence on the viewing position is small and can be neglected.

3. As has been mentioned in the previous section, the irrad~:nce EeS of each point (x,y) in the image plane R 2 at a viewpoint teR is assumed to be roughly proportional to the radiance of the point in the scene projecting into the point in the image plane at the viewpoint t, and the proportionality coefficient, to be independent of the location and the viewpoint t.

4. The scene is assumed to be made out of opaque objects. The criteria for selecting feature points and the criteria for assigning feature values to each feature point M are I ~I 28 66,645 assumed to represent the true properties of the objects and to be independent of the objects' spatial attitudes.

The velocities of neighboring points of the objects in the scene are assumed to be similar, on one hand, and to change slowly with respect to viewpoint, on the other hand. In other words, it is assumed that the parametric viewpoint-varying velocity vector field of the surface points of each object in the scene varies smoothly across the space and across the viewpoint location.

Let lt be a subset of the space R 2 defined by the relation t j(A t,A t )eR 2 ,A t+A t 0, (t+At E" tf i

K

I

(3-1) let Pt be a subset of the set It, and let G t G x Pt be the 15 Cartesian product of the set G on the set Pt.

The elements from the set G t are denoted by the symbol gt, and the set G t is assumed to be a measurable space with some measure dg t defined on the family B(Gt) of Borel subsets of the set G t In practice, the set G t is finite, and the 20 measure dg t is a point measure on the finite set Gt, but a more general approach is used for the sake of uniformity of A.s presentation.

Let gt (g,At',At*)eGt,aC[l,o) and let gtO(x,y,t,u',v a be a function defined by the relation g, y, t, u, v (g (X+A t y, t) ,y+A t t+A t) -go (x-A t-uO(x, y, t) ,y-A t- t-A (A t-+A for every (x,y)en, where: is the velocity vector of the image corresponding to the parameter value a at the image point and the viewpoint t. Here, to simplify the notation, the symbol gt is used as the part of the notation gt for the above defined function in addition to its role as the index for such a function.

I rd I~ li~YP~C~ -a~ll- 29 66,645 Let (x,y)en be an image point that is taken at the viewpoint t projection of some point in the scene that does not belong to the occluding boundaries. The occluding boundary of an object is defined as the points in the scene belonging to the portion of the object that projects to its silhouette. The assumptions 1-4 made at the beginning of this section imply that the absolute value of the function gtJ(x,y,t,uO,vO) is small. Therefore it is natural to use the function g? y, t,1u, v 0 2 (3-3) as a part of the functional whose minimum specifies the estimate a, 0 of the velocity vector of the image corresponding to the 15 parameter value a. It is customary for the function (3-3) to be called the "optical flow constraint" corresponding to the parameter value a and the index gt 12-14, 18-21].

Note that if the viewpoint increments At', At approach zero, then for an integer n greater than or equal to 2 the n th-order partial derivative of the function gto(x,y,t, O 0O with respect to the components ao a of the estimate of the velocity vector approaches zero at a rate proportional to n 1 This, in particular, implies that when the viewpoint increments At',At approach zero the function gt 0 (x,y,t,1, o approaches the function which is linear with respect to the components uo,I of the estimate of the velocity vector.

It can be easily verified that the estimate of the velocity vector on which the function achieves its minimum satisfies the system of equations go(x,.y o O )g ix,y, t, o gi y, t, gI t, 0 0) =0, (3-4) where the functions gto o 7o0 and gta o o V o are i 4 P l ill C- -PII( 30 66,645 the first-order partial derivatives of the function gto(x,y,t,a,9 a with respect to the components i and -V of the estimate of the velocity vector.

Let the vector (Au 0 be defined as in the relations z2 y, t) =u a 0) +Au y, V° +Av° t).

(3-6) By expanding the left hand sides of the relation into 0 their Taylor series about the velocity vector and by setting the second-order and higher partial derivatives of the function gt"(x,y,t,u,v with respect to the components u° ua(x,y,t) and v° v(t) of the velocity vector to zero, we obtain the system of linear equations, relative to the vector t(glu (xy, t,u v 2A t) +g t, 0

O)

g(X, y, t, u, vO A +gu t, u a v a g y, t, u, v) =0, gFu(xy, y, tU V g Y (XY, y, U, uy) A uo(x, y, t) g v =0.

eo* (3-7) where the functions gtu°(x,y,t,u,v a and gtv (x,y,t,u ,v 0 are the first-order partial derivatives of the function gto(x,y,t,u,vO) with respect to the components u° and vO of the velocity vector.

Let us first assume that the norm of the gradient -ap re -P III~ CR~ ~sB~ II I- 31 66,645 (9gu y, t, u vO) gt y, t u v)) (3-8) of the function gtO(x,y,t,uI,v) with respect to the components u" and v a of the velocity vector, is relatively large; then the system of equations is equivalent to the relation gtu y t, UC, Vu Au (x y, t) +gtv(x,yft, U t g t, u, v 0 (3-9) 0 The last relation implies that the requirement for the function to achieve a minimum on the estimate of the velocity vector o specifies the component of the estimate in the direction of the gradient and leaves the component orthogonal to the gradient undetermined.

If the norm of the gradient is relatively small; then the system of equations becomes too weak, and the requirement for the function to •achieve a minimum on the estimate (ao(x,y,t),Q 0 of the velocity vector imposes no constraint on 20 the estimate.

The above discussion suggests that we need to impose additional constraints on the estimate 0 (x,y)en, ae[l,o) of the parametric velocity vector field (x,y)En, e[l1,) in order for the computation to be well defined.

Let S=((sx,Sy) I sy) eR2,x2+Sy2=1) be a unit circle in the two-dimensional Euclidean space R 2 with each vector sm(s,sy)ES specifying the direction in the image plane. The set S is assumed to be a measurable space with some uniform measure ds defined on the family B(S) of Borel subsets of the set S.

Let s=(sx,,y)eS, and let (s,VuO(x,y,t)) be the function defined by the relation i. II 32 66,645 (s,Vu a t) x+sx,y+Osy, t) O, (3-10) for every (x,y)en. Here the symbol V stands for the gradient of a function with respect to the variables x,y.

Let s (sx,Sy)ES, and let (x,y)en be an image point that is taken at a viewpoint t projection of some point in the scene that does not cross the occluding boundary in the direction s, meaning that for all sufficiently small positive constants o, the following two points in the scene both belong to the same object and are sufficiently close to each other: the first is the one projecting to the image point (x+-os,y+Osy) at the viewpoint t, whereas the second is the one projecting to the image point at the viewpoint t. The assumption 5 made at the beginning of this section implies that the absolute 15 values of the function is small. Therefore it is natural to use the function o V y, t)) 2 (3-11) :as a part of the functional whose minimum specifies the S 20 estimate of the velocity vector of the image corresponding to the parameter value a. It is customary for the function (3-11) to be called the "smoothness constraint" in the direction s[2,12-14,18-21].

Finally, to make the process of computing the estimate 0 of the parametric velocity vector field (x,y)en,ae[l,oo)) well defined, we require that for every image point (x,y)en and for every parameter value ae[l,o) where the optical flow constraints and the directional smoothness constraints impose no restrictions the difference between the unknown estimate 0 of the velocity vector (ua(x,y,t),vo(t)) and its initial estimate (0oo(x,y,t),90a(t)) be small. This can be accomplished by including the functions IP~l~l~ I 33 66,645 (u (Xy, t) -I y,t))2, (3-12) -9 2 (3-13) in the functional whose minimum specifies the estimate of the velocity vector of the image corresponding to the parameter value a in addition to the functions and (3-11) specifying the optical flow and the directional smoothness constraints, provided that the weights associated with the functions (3-12) and (3-13) are small relative to the weights associated with the functions and The initial estimate (uo 0

,Q

0 appearing in the relations (3-13) is defined later in this section. The constraints imposed 15 on the estimate of the velocity vector at the image point (x,y)en and the viewpoint t by the above requirements will be called the "regularization constraints".

As has been stated earlier, the optical flow 20 and the directional smoothness constraints are not necessarily valid at the points near the occluding boundaries, even when the assumptions described at the 9 *beginning of this section are observed. The method of computing the estimate of the parametric velocity vector "o*o 25 field of the present invention resolves the above difficulties by adjusting the weight associated with each constraint in such a way that it becomes small whenever the constraint is not valid. In the following, it is described how the functions and specifying the optical flow, directional smoothness, and regularization constraints, are combined into the functional of the estimate of the parametric velocity vector field. The estimate is then computed by solving the system of nonlinear equations arising from a certain optimality criterion related to such functional.

The following describes the system of I r rrcv la 34 66,645 nonlinear equations used in the method of the present invention in its estimation of velocity vector fields. In the method of the invention, the estimate a I (x,y)esnfl, ae[l, of the parametric velocity vector field (uO(x,y,t) I Enfl,ae[l,)) is determined as the parametric vector field on which a weighted average of the optical flow, directional smoothness, and regularization constraints is minimized.

The weight functions are chosen in such a way that at every image point (x,y)E which is near the occluding boundary the following requirements are met: each weight function associated with an optical flow constraint becomes small whenever the optical flow constraint is not satisfied, and each weight function associated with a smoothness constraint corresponding to a direction in which an occluding boundary is crossed becomes small whenever the directional smoothness e.

constraint is not satisfied. The criteria for the presence of the occluding boundary near a given image point, and the criteria for the occluding boundary being crossed near a 20 given image point in a given direction, can be defined most effectively in terms of the values of the unknown estimate of the parametric velocity vector field. Therefore, the values of the unknown estimate of the parametric velocity vector field have to be used implicitly in each of the S. 25 weight functions. On the other hand, each weight function has to be treated as if it were independent of the values of the unknown estimate of the parametric velocity vector field, because it only specifies a relative significance of the corresponding constraint as a part of the weighted average and not the constraint itself. To overcome these difficulties, two copies of the unknown estimate of the parametric velocity vector field are introduced: the invariable one, and the variable one. In the weight functions the values of the invariable copy of the unknown estimate of the parametric velocity vector field are used, whereas in the constraint functions the values of the variable copy of the unknown estimate of the parametric velocity vector field are used. Then a more general variational principle as opposed to an energy minimization i IY~I__ 35 66,645 principle is applied to the weighted average of the optical flow, directional smoothness, and regularization constraints to derive the system of nonlinear equations of the unknown estimate of the parametric velocity vector field. This variational principle can be described as follows.

Let (d,V) (Ux,y,t) 0 I (x,y)en,ae[l,oo)) be an invariable copy of the estimate of the paramet: ia velocity vector field, let (x,y)En,ae[l,oo)} be a variable copy of the estimate of the parametric velocity vector field, and let the parametric vector field (AO,Av) be defined by the relation (A ,A -(u 7 (x c e 1, (4-1) Assume that the derivative of the parametric i vector field (AO,Av) in the direction normal to the boundary an of the image n is equal to zero for every parameter value ae[l,oo) and for every image point belonging to the 20 image boundary an. This will impose restrictions on the types of the variations that can be applied to the variable copy of the unknown estimate of the parametric velocity vector field.

Let be a functional of the 25 parametric vector fields defined as a weighted average of functions and (3-13), specifying the optical flow, directional smoothness, and regularization constraints, respectively, by the following relation: 9 1- s8 -36 66,645 I. a Gt: t) t) (4-2) Here: aci)'Y xyt' rloV I 60,ae[(1, 10tGt is a C weight associated with an optical flow constraint corresponding to the parameter value a and the index gt, which is a function of the independent variables x,y,t and *of the f unctions y t) 0 a 0 v ai oI IvQ(X, y, t) I of these variables; PS o(X, y, t, jYaa, Vf 0 (x,y)efl,ce[1,oO),scS is a weight associated with a smoothness constraint in the direction s, which is a function of the independent variables x,y,t and of the functions Eua~xy~t)vO~t,(suO) (s,VQ of these variables; 15 (ii y, 0 are positive constants specifying weights for the regularization constraints; (iv) the function IVQ'(x,y,t) I ,(x,y)efl,ae~l,co) is the norm of the vector valued function Vil 0 which is the gradient with respect to the variables x, y of the function fi1(x,y,t), given by the relation (4-3) the functions (x,y)efl,a6(1,co),seS are 37 66,645 defined as in the relation (3-10); (vi) the function gta(x,y,t, (x,y)en,ae[l,o) ,gteGt is defined as in the relation and (vii) the components o,1 0 o0(t) of the initial estimate of the velocity vector corresponding to the parameter value oa[l,o) at the image point (x,y)en are defined later.

The estimate of the parametric velocity vector field is then defined as the parametric vector field on which the functional considered as the function of the parametric vector field and depending Son the parametric vector field as on the parameters, achieves a local minimum when the value of the parametric vector field is identically equal to the value of the parametric vector field Taking into account the "7 relation the functional f(i,v,u,9) can be expressed in the form The parametric vector field S' (AO,Av) specifies a perturbation to the parametric vector 20 field and the functional f(u,9,u+Au,-+Av) assigns the cost to each choice of the parametric vector field and its perturbation (Af,AV). Then the estimate of the parametric velocity vector field is the parametric vector field for which a locally minimal cost is achieved 25 when the perturbation is identically equal to zero.

Using the calculus of variations with respect to the parametric vector field (Af,,AV) applied to the functional f(u,V,O+AO,v+Av), the above defined estimate of the parametric velocity vector field is a solution of the system of equations ar 38 66,645

S

(4-4) fffc (X,y ,Li 0 1 Iva gtO(X, Y, t,U0i t)g where the functions gt~a(X,yt,a,,rP) and gtO 0 (x,y,t,Qe 1 'Qe) are the first-order partial derivatives of the function gtO(x,y,t,f1O,,V) with respect to the components fi i(x,y,t) *and r iz(t) of the estimate of the velocity vector.

To select the weight functions associated with the optical flow, directional smoothness, and regularization constraints, consider the functional fi,'i,)defined by the relation undor the conditions that the parametric vector field (i0,4) be identically equal to the parametric vector field and Utia both be equal to the estimate (QO y,t) ,V y) en ac C1, z) of the parametric velocity vector field obtained as the solution of the system of nonlinear equations In this case the functimnad f( takes the f orm 39 66,645 f (11,VIl) =f (ff fO0.5r (XiY, t,Q a, V,G, jVOa I) 1 0 Gt, VOY(X, y, t) 2 dS+ 0 UG(t) Vol' 2] doa.

Let ((lo y, t) "00(t))j(x Iy) cl, a [Ioo) be an estimate Ct' the parametric velocity vector fk.,which is obtained as the solution of the system, of nonlinear equations Select each of the weiciht functions ago(,yttIa 0"11VgI),qi~ in such a way that the contribution of the optical flow constraint gt 0 to the functional becomes small for every image point cO located near the occl,>,ding boundary where the :optical flow constraint gt 0 (xy,t,Tl 0 1 is not satisfied.

Similarly, each of the weight functionz p 8 eS will be chosen so that the contributions of the directional smoothness constraint (s,VfIO(x,y,t)) to the functional become small for every image point (x,y)en l ocated near the occluding boundary where the occluding boundary ii5 crossed in the direction s and the directional smoothness constraint (s,VCIY(x,y,t)) is rnot satisfied. Let the image point (x,y)cnf be near the occluding boundary; then the following' two ver~its are likely to happen: 1. For some small positive viewpoint increments At' and At" the point in the scene projjecting into the image point at the viewpoint t either is visible at the viewpoint and is invisible at the II C_ 40 66,645 viewpoint or it is visible at the viewpoint (t+At and invisible at the viewpoint Let (xl,y 1 and (x 2 2 be projections of such a point in the scene into the image plane taken at the viewpoints (t-At) and (t+At respectively. If the radiance of the point in the scene projecting into the image point (xl,yl) at the viewpoint is significantly different from the radiance of the point in the scene projecting into the image point (x 2 ,y 2 at the viewpoint and if the radiance undergoes some changes in the neighborhoods of the above points in the scene, then the following cases are likely to be encountered: 1.1. The point in the scene projecting into the image point (x-At'f(x,y,t) ,y-At' 0 at the viewpoint and the point in the scene projecting into the 'image point at the viewpoint both belong to the occluded object, on o-l hand, and have different radiances, on the other hand The later may happen if the above points occupy distinct locations on 20 the object of the scene which, in turn, may be caused by the presence of the occluding boundary. Under these conditions the absolute values of some of the functions gOt(x,y,t, o a,

O

,g t eG t become large.

1.2. The point in the scene projecting into 25 the image point (x At'f(x,y,t),y-At-f(t)) at the viewpoint and the point in the scene projecting into the image point at the viewpoint both belong either to the occluding object or to the occluded one. In addition to this, for some not-too-large vector (Ad°(x,y,t),Av90(t)) one of the following two points in the scene belongs to the occluding object, while the other belongs to the occluded one: the first point is the one projecting into the image point (x-At'( o (x,y,t)+AaO(x,y,t)),y-At (9C(t)+Ao 0 at the viewpoint and the second point is the one projecting into the image point (x+At+ o (x,y,t)+Af o O (t)+AV o at the viewpoint Likewise, for the same vector one of the following two points in the b -~IC~L I I Y_ 41 66,645 scene belongs to the occluding object, while the other belongs to the occluded one: the first point is the one projecting into the image point (x-At(fa(x,y,t)-AIO(x,y,t)),y-At at the viewpoint and the second point is the one projecting into the image point (x+At (oa(x,y,t)-Af at the viewpoint (t+At) (see Figure In such a case the values of the functions (gt o (x,y,t,f, o g) 2 ,gteGt become relatively small, while the values of some of the functions (gto(x,y,t, fa+Aa, a+A))2',gtGt, as well as the values of some of the functions (gta(x,y,t, la-Aloga-AC o 2, gteGt, become O relatively large. This implies that the values of some of the functions 1* (g t, a t +Aa o 1 0 +Af) 2 tao-AA o Q, a -A

O

2 -2 g y, t' ila, T) 2 (4-6) gteGt become large. If the second-order and higher partial derivatives of the function gt o (x,y,t, 0 f 1 with respect to the components a=la(x,y,t) and of the estimate of the velocity vector are ignored, each a£ the functions (4-7) can be expressed in the form 2 (A y, t) go t, D, 0) +A (t) g t, a° 2 (4-7) where the functions gtoa(x,y,t, aa, a) and gt o d, a) are the first-order partial derivatives of the function gt"(x,y,t,Ca',o) with respect to the components ii and tO of.

the estimate of the velocity vector. For the vector s s (sx,sy), seS that is collinear to the vector (Afo(x,y,t),Ava(t)), the above observations imply that the absolute values of some of the functions (s,V'gto(x,y,t, C, 0 defined by the relation ,I IL~a~b 42 66,645 Vg y, tC =Sxg (xy, t, O, vi) sygto(xy, t, u, 0 (4-8) become large.

2. The estimate ((ft(xy,t) (tx,y)efn,ac[(l,co) of the parametric velocity vector field changes rapidly in the vicinity of the image point (see Figure which implies the following: 1 e 2.1. The values of the functions jVIa(x,y,t) become large.

2.2. The absolute values of the functions (s,VV(x,y,t)),seS become large, provided that the image point crosses the occluding boundary in the direction s.

S.

15 If the conditions on the scene being imaged and on the imaging process itself described at the beginning of the previous section are satisfied, then in the event of the absolute value of the function gto(x,y,t,aO,9V) being S* large for some gteGt, or of the values of the functions IVQ(x,y,t) being large in addition to the absolute value of the function gt o a being large for some geGt, the image point is likely to be near the occluding boundary, while in the event of the absolute values of the •functions (s,Vda(x,y,t)),seS being large, or of the absolute values of some of the functions (s,V'gta(x,y,t,o, 0 0 gteG being large in addition to the absolute values of the functions (s,VQ 0 seS being large, the image point is likely to cross the occluding boundary in the direction s. Note that in the case of the functions |VQ(x,y,t) j being large, the image point does not necessarily lie near the occluding boundary. It may, for example, lie on the object of the scene whose projection onto the image plane undergoes a rotation or a local deformation. Also note that in the case of some of the functions (s,V'gt(x,y,t, gtEGt being large, the image i id I~IIISs lP""I~I~LIPI~I 43 66,645 point does not necessarily cross the occluding boundary in the direction s. It may, for example, cross the radiance boundary arising from a texture or a sudden change in the illumination.

These observations suggest the following: each of the weight functions a0(x,y, t, Q, IVuI1j),gEG should be a steadily decreasing function relative to the absolute value of the function gt o (x,y,t,Cao,Q') and to the values of the function IVfo(x,y,t) multiplied by the absolute value of the function gt

O

(x,y,t,uo,Qo) and each of the weight functions p'o(x,y,t,lao,, (s,V a seS should be a steadily decreasing function relative to the absolute values of the function (s,VaO(x,y,t)) and to the absolute values of some of the functions V'gta(x,y,t,o, eGcG multiplied by 15 the absolute values of the function It is natural to require each of the weight functions corresponding to either an optical flow constraint or a directional smoothness constraint to be selected in such a way that the product of the weight function and the corresponding constraint function appearing as the summand Sin the functional becomes larger whenever the constraint function becomes larger. This makes the solution of the system of nonlinear equations determining the estimate of the parametric velocity vector field more robust with respect to the choice of the initial estimate of the parametric velocity vector field and with respect to the noise. To be more specific, these requirements consist of the following.

1. For every gteGt the function defined by the relation: at(X 'y tx l, a, |VO °|ll (x y, Q, o o 2 t (4-9) and appearing as the summand in the functional is required to be steadily increasing relative to the function (gt o 2 In other words, the function is PIS~PS 1II~RPII Cls 44 66,645 required to increase whenever the function (gto(x,y,ti,.Y,'Q)) 2 is increased.

2. For every seS the function defined by the relation: (4-10) and appearing as the summnand in the functional is required to be steadily increasing relative to the function (s,Vi3L(x,y,t)) 2 In other words, the function (4-10) is required to increase whenever the function 2 AO is increased.

W The weight functions 9 (4-11) 15 for every gt~t where r 2 p 2 q 2 are non-negative real constants, and (a I+ (C 2 +b I(x, V gg(x, y, t fj 2 (S Va y 2 -1 (4-12) for every seS, where the expression is def ined as in the relation (4-13) and a 2 C 2fb ;12 G ar no-negative real constants, comply with the above requirements and with the requirements stated _IC 45 66,645 earlier in this section.

The approach of obtaining the estimate of the parametric velocity vector field (x,y,t,Q 0 (x,y)en,ae[l,oo) by solving the system of nonlinear equations for each parameter value a independently may encounter difficulty relating to the fact that the system of equations may, in general, have many solutions. To overcome this difficulty, the values of the estimates of the velocity vector field corresponding to the different values of the parameter a are tied together by imposing the following restriction: the estimate of the parametric velocity vector field S{ I en,ae[l,o) is required to be a continuously differentiable function with respect to the parameter a for every image point (x,y)ef and for every parameter value ae[l,o). The additional restrictions on the estimate of the parametric o. velocity vector field 0 j en,aE[l,) are imposed in the form of the boundary conditions as 20 follows: 1. The estimate of the parametric velocity vector field 0 I en,ae is required to converge to the vector field that is identically equal to a given vector constant when the parameter a converges to o, and such convergence is required to be uniform over the set of the image points n.

2. The derivative of the estimate of the parametric velocity vector field (x,y)en,ae[l,o)) in the direction normal to the boundary an of the image n is required to be equal to zero for every image point belonging to the image boundary an and for every parameter value ae[l,c).

Let ae[l,o) be a given parameter value, and let Aa be an infinitesimal positive real value. Taking into consideration the above restrictions, the initial estimate of the velocity vector field {(0o 0 (x,y,t),9 0 (x,y)En) corresponding to the parameter value a can be defined as follows:

I

I~IB~BP~~ 1 46 66,645 a (Xy, t) 9 Q* y, y)eQ (4-14) Then the regularization constraints (3-12) and (3-13) can, respectively, be approximated by the relations Yu,o y, t) -IOU y, o y (4-15) Yv,o 0 t) i a (4-16) for every image point (x,y)Enf, where YuYv are positive real constants given by the relations Yu=Yu,oAG.

(4-17) Yv=Yv,o YV"=y 0 Aa.

(4-18) 15 By substituting the regularization constraints Yu,0( C(x,y,t) and Yv, 0 appearing in the system of nonlinear equations with the constraints a a -Yua- 0 t) and -y o respectively, and by aa *aa defining the weight functions a" t, I VO °),g;ce and P/3(x,y,t,u°,Q 0 according to the relations (4-11) and respectively, the system becomes the following system of equations relative to the unknown estimate of the parametric velocity vector field I en, ae[ C _i I ~-161 L~Y 47 66,645 G, r+ (p q1 VCla1(x, y, t) I2) Y, tU 0 19'ry a 0 f V 2 V71 a(.x)Y=0 d=.

constraints b (s5y, y, t, 0 0j,y 2 an5~( 0 t ff 0 2 pprn 22In[ 1(,y g(,y the fucioa 0~i,4f/O) 2ddefndy the elati nd-2 ereplacedion th the constrais contrnt 0.5y(Q(x yu,ci(x, y,t) flo yt) 2 and 0 y~ V t 0. 5yv(C-(t) 2 respectively, and provided au the weight functions ao (XOYl t:1FC7 a,,Vacy j gG. are defined according to the relation whereas the weight functions paxyoa (s,V 0 SES are defined according to the relation the functional of the parametric vector field (ii,Q) (la(xyt) Cr~t y) e n,a c[1, o) I and the parametric vector f ield x O t, I(,y ea)[1 o defined by the relation can be replaced with the functional h(a~~fQ of the same parametric vector fields 48 66,645 defined as in the relation h(g, y, 2) J V if ln(r 2 (p 2 q

V

a x ,y, t 12 (gV(xy t -dg))2) 1 V 2 (p 2 +q 2 V[z j t) 12) fin(a2+(c 2 +b 2 (sV g I(x,y, t, a, V1Q(x,y, t)) 2 t ds S2 (c 2 +b 2 Vg t j, O) 2) t) t) a (Xy t) 2] dxdy+ aa 0 5 Y -e 0 1 t) f do, Ldo, 00 (4-20) 'where the symbol in stands for the natural logarithm 5 function.

*ee* The above substitutions imply that the •parametric vector field (fi,v) satisfies the system of •nonlinear equations (4-19) if and only if every variation of the functional h(u with respect to the parametric 10 vector field is equal to zero whenever the parametric vector field is identically equal to the parametric vector field P' Note that when the constant q 2 approaches zero and when the increments between successive viewpoints of the 0**0 15 viewpoint sequence uniformly approach zero the functional approaches the functional, which is independent of the parametric vector field and the problem of finding the estimate of the parametric velocity vector field becomes the energy minimization problem.

Let en,a£[lo)) be the estimate of the parametric velocity vector field obtained as the solution of the system of nonlinear equations We now examine the behavior of the system of nonlinear equations (4-19) near such a solution. First, the changes in the behavior of 49 66,645 the optical f low and the directional smoothness constraints when the image point goes from the position beyond the occluding boundary to the position near the occluding boundary are considered. Then, the dependence of the optical flow and the directional smoothness constraints on the parameters r 2 p 2 q 2 specifying the weight functions ag o y't' ia,CaJIV1a) gjGt as in the relation (4-11) and on th praetrsa2' C22, gtFea specifying the weight functions '8'a t' a(S, vfia) seS as in the relation (4-12) are investigated.

Let (AiO,A'O) {(Afi" y,t) A~r(t) I y)Efl, cc[l1, be a parametric vector field defined by the relation then in the neighborhood of the parametric vector field the functional can be expanded into the quadratic form of the parametric vector field as follows: (rff (p 2 +q 2 jVjO t) 12) 2 (gt? 0 (x'-yt riO 0 1 fO) AU'(x'Y, 0 (C 2 +b 2 V'g I y, t, fj 2) V, 1g, y, 2) 2 V(A y' t) 2 ds 4y(A tOa(X,y_, t) 2 dXdy+y, t) 2 dca.

(4-21) IISBa~;meaPlrrrssPlsIrP-ls-~ 50 66,645 The optical flow constraint gta(x,y,t,f',9 0 corresponding to the index gtEGt is represented in the system of nonlinear equations (4-19) by the expression (g t, 0 2 r 2 +(p 2 +q 2 jVj"(xy, t) 2 (gg (xy, t, Q 2 (4-22) With respect to this expression the following cases can be considered: 1. The image point is beyond the S occluding boundary. Then the value of the function (gt o (x,y,t,f o ,9i)) 2 may become small. In this case the value of the function (4-22) would approach the value of the function r e* S(4-23) This, in turn, would make each, of the two summands of the system of nonlinear equations (4-19) associated with the opt'cal flow constraint gto(x,y,t, i corresponding to the index gt behave as a linear function relative to the function grt(x,y,t,Q,CO) with the proportionality coefficient equal to the function gto(x,y,t,QO,V)/r 2 in the case of the first summand and to the function gt c(x,y,t,~ac,Q)\r 2 in the cas- if the second summand.

2. The image ,int approaches the occluding boundary. Then the value of the function (gta(x,y,t,o 2 may increase. In the case when the value of the function (gt(x,y,t, a,o)) 2 is less than the value of the function r 2 /p 2 q 2 jVI o(X,yt) 12), the fact that the summand of the quadratic form (4-21) associated with the optical flow constraint gtO(x,y,t,fi,o) is positive implies that the function (4-22) would increase. In the case when the value of the function (gt o (x,y,t,fo,Qo)) 2 is greater than i _I rr 51 66,645 the value of the function r 2 /(p 2 q2JVo(x,y,t) 12), the fact that the summand of the quadratic form (4-21) associated with the optical flow constraint gt (x,y,t,O,'O 0 is negative implies that the function (4-22) would decrease. Therefore, the function (4-22) achieves its maximal value 1 2r/p 2 -+q 2 (V y, t) 2 (4-24) when the value of the function (gt o 2 is equal to the value of the function r 2 /(p 2 q21VV(x,y,t) 12).

It may also happen that the values of the function (Vo(x,y,t) becomes large. Under these conditions the following values are reduced: the maximal value of the function (4-22) specified by the relation the value r 2 /(p 2 q21VrV(x,y,t)12) of the function (gto(x,y,t,1o,Va 0 2 where this maximal value is achieved. The above situation may occur in the following cases.

1 The occluding object may artificially constrain the estimate 0 in the direction orthogonal to the occluding boundary. Then the optical flow 20 c nstraint may react by shifting the estimate of the velocity vector in the direction parallel to the occluding boundary to reduce the absolute value of the function gt o (x,y,t,a,V o This, in turn, may increase the value of the function JVa(x,y,t)|.

25 2. The estimate may diverge as the result of nonlinearity of the cptlcal flow constraint or as the result of the noise present in the images.

The directional smoothness constraint (s,V1a(x,y,t)) in the direction seS is represented in the system of nonlinear equations (4-19) by the expression (s,Vg(x,y, t))2 a 2 +(c 2 +b (s,Vgg(lx,y, to,V-o 2 (s,VQo(x,y, t) 2 (4-25) With respect to this expression the following cases can be II N il~P-- r 52 66,645 considered: 1. The image point does not cross the occluding boundary in the direction s. Then the value of the function (s,VY(x,y,t)) 2 may become small. In this case the value of the function (4-25) would approach the value of the function s, Vg° t)) 2 a 2 (4-26) This, in turn, would make the summand of the system of nonlinear equations (4-19) associated with the directional smoothness constraint behave as a linear function of the corresponding constraint with proportionality coefficient equal to the cor.tant value 1/a 2 15 2. The image point approaches the occluding boundary in the direction s. Then the value of the function (s,VO(x,y,t)) 2 may increase. In the case when the value of the function (s,V l 2 is less than the value of the function a2/(c+b2(sVgtO(xt, 0 the fact that the summand of the quadratic form (4-21) associated with the directional smo( Thness constraint is positive implies that the function (4-25) would increase. In the case when the value of the function (s,Vua(x,y,t)) 2 is greater than the value of the function (C -b 2 (s,V'gto(x,y,t,avo))2), the fact that the summand of the quadratic form (4-21) associated with the directional smoothness constraint is negative implies that the function (4-25) would decrease. Therefore, the function (4-25) achieves its maximal value 1 (4-27) 2ac+b2 V t, )2 when the value of the function 2 is equal to the value of the function a 2 /(c 2 +b 2 (s,V gto(x,y,t,'aga))2).

In the case when the values of some of the I i -II rslPs~-Y L~l~sl L~Cr~P-~ P- I 53 66,645 functions (s,V'gt o (x,y,t,l",gr)),gteG are 7.arge the following values are reduced: the maximal value of the function (4-25) specified by the relation the value a 2 /(c 2 +b 2 (s,V'gto(x,y,t, o ga 2 of the function (s,Voa(x,y,t)) 2 where this maximal value is achieved. Such an action reflects the fact that in this case the image point is likely to cross the occluding boundary in the direction s if the value of the function (s,Vi7(x,y,t)) 2 specifying the degree of the smoothnesm of the estimate of the velocity ve-tor field in such a direction s is large.

has been stated earlier the role of the regularization constraint is to discriminate between the 9 parametric vector fields giving the same optimal values to the weighted average of the optical flow and the directional smoothness constraints without causing any significant changes in these values. This is achieved by assigning small values to the parameters Yu and y, appearing in the system of nonlinear equations For the sake of simplicity of the analysis given below, we shall ignore the regulari ation constraints by assuming that the parameters Yu and y, are equal to zero while keeping in mind that the solution of the system of nonlinear equations (4-19) is locally unique.

Taking the above into account, the process of solving the system of nonlinear equations (4-19) described in the next section can be regarded as an attempt to achieve an equal balance between the functions representing two antagonistic sets of constraints: optical flow constraints :and directional smoothness constraints. This, in part, is the consequence of the fact that the quadratic form arising from the linearization (5-16) of the system of nonlinear equations (4-19) described in the next section is positive definite.

It is not difficult to see that by increasing the values of the parameters r 2 p 2 q2 apnearing in the relation we decrease the values of the weight functions a y, C, JVoi) gf.Gc This would shift IIE~Ps~8~P'--a 54 66,645 the balance in the system of nonlinear equations (4-19) toward the directional smoothness constraints. To restore the equal balance between the optical flow and the directional smoothness constraints the system of nonlinear equations (4-19) would react in such a way that the absolute values of the function (s,V~ 0 (x,y,t)),seS corresponding to the directional smoothness constraint is reduced, while the absolute values of the functions gt(x,y,t, 0 ,I ),g t eG t corresponding to the optical flow constraints are increased.

Likewise, by increasing the values of the parameters a 2 c 2 b 2 ,geG, appearing in the relation we decrease the values of the weight functions

P

8 s(x,y,t,°,v 0 (s,Vl)),seS. This, in turn, would shift the balance in the system of nonlinear equations (4-19) toward the optical flow constraints. To restore equal balance between the optical flow and the directional smoothness constraints the system of nonlinear equations (4-19) would react in the way that would reduce the absolute values of the functions gtO(x,y,t,,O) g t e G t corresponding to the optical flow constraints and would increase the absolute value of the function (s,VUt(x,y,t)),seS corresponding to the directional smoothness constraints.

Taking into account the earlier discussion related to the behavior of the optical flow and the directional smoothness constraints, we can conclude the following. The parameters r 2 and a 2 have their greatest impact on the solution of the system of nonlinear equations (4-19) at the image points that are away from the occluding boundary, while the parameters p2,q2 c2,b ,ge, geG mostly influence the solution at the image points that are near the occluding boundary. The parameter r 2 defines the coefficient of proportionality for the functions gtO(x,y,t,"4,9 0 ),gt c G t specifying the optical flow constraints, while the parameter a 2 determines the proportionality coefficient for the function s 55 66,645 specifying the directional smoothness constraints. The combination (4-24) of the parameters r 2 p 2 and q 2 determines the upper limit for the influence of the optical flow constraints on the solution of the system of nonlinear equations while the combination (4-27) of the parameters a 2 c 2 and b 2 gG, determines the upper limit for the influence of the directional smoothness constraints on the solution of the system of nonlinear equations (4-19).

Earlier the estimate 0 I (x,y)en,ae[l,o)) of the parametric W velocity vector field I en,cre l,o) was defined as a solution of the system of nonlinear equations which can be expressed in the following form: G r 2 +(p 2 2 a02) 0 2 f(s, V Vi) ds s a 2 +(c2+b Vg)2) Vua) 2 f: -y dgodxdy a-Q r 2 +(p 2 +q 2 V 2 2) )2 (5-1) Here, for the sake of compactness of the notation, the arguments of the functions appearing in the system of nonlinear equations (4-19) have been omitted. Also, the symbols gt0, §tg 0 4t have been used instead of the symbols gt°, gto, gt° to indicate that at every image point (x,y)en the corresponding functions depend on the estimate of the velocity vector i Ic-aF" ~p~a~ill~a~aar~8 56 66,645 To compute the estimate of the parametric velocity vector field 0 I a restriction is imposed ourselves to a finite subset of the parameter values from the interval Let a 0 be the parameter with the value equal to o. Suppose that the value of the parameter a, is equal to 1; and that the sequence of the parameters ao, a l n is not too sparse, and is decreasing in values. Then the estimate of the parametric velocity vector field j a n is obtained by successively solving the system of nonlinear equations (5-1) for each value of the parameter a=a 1 an starting with P the parameter value a=ao where the solution is identically equal to a given vector constant and proceeding in the direction decreasing with respect to the parameter values using the solution of the system of nonlinear equations obtained for the preceding parameter value as the initial estimate.

Let 6(a) be a function of the parameter a 20 defined on the set of the parameter values a=al,. ,a by the following relation: (5-2) If the weight constants yu,yv, 0 are defined as in the relations Yu Yu( (5-3) Yv (5-4) then the regularization constraints -Yuo-- and a -YV aI f can respectively, be approximated as -I I III sr~C-_~LsC~ 57 66,645 t) y, t) y, t)) -Y y V (V t) -'VI a t, ao (5-6) for every image point (x,y)en, and for every parameter value a n Under these conditions the system of nonlinear equations of the unknown estimate of the parametric velocity vector field en,a=a, ,an Stakes the form f r2+(p 2 +q 2 V a 2 (g)2 S0 (s,V Va) ds+ 1 0 s a 2 +(c 2 +b (s,VIgg) 2 (s,Vu°) 2 yO (a _0 (oa) =0, r2+ (p2+q2 Vi 2) 2 (5-7) Then the computational process for estimating the parametric velocity vector field can be more precisely described as follows: for every parameter value a=a, using the previously obtained estimate ,9 6 I en) of the velocity vector field corresponding to the parameter value 6(a) as the initial estimate, find the estimate I (x,y)en) of the velocity vector field corresponding to the parameter value a that is the closest to the initial estimate by solving the system of nonlinear equations For a fixed value of the parameter sl -1 3 Is_ ~e ra ~F-n~o 58 66,645 a=a 1 the system of nonlinear equations (5-7) relative to the unknown estimate (x,y)En) of the velocity vector field corresponding to the parameter value a, defines, for every image point (x,y)en, a set of relations between the components of the velocity vector and its partial derivatives, which can be expressed in the form

F°(U

0 0.

(5-8) To find the estimate (Ol,g

C

of the velocity vector field corresponding to the parameter value a as an estimate of the solution of the system of nonlinear equations (5-7) corresponding to the parameter value a that is the closest to the initial estimate 5 the iterative updating scheme can be used. In this scheme, starting with the initial estimate (e6(a), 6 which is taken as the current estimate of the solution of the system of nonlinear equations corresponding to the parameter value a, the improved estimate 0 j n} of 20 the solution of the system of nonlinear equations (5-7) corresponding to the parameter value a is determined by the step (u a 8(A°,AV), (5-9) where the scalar parameter defines the length of the step, while the vector field (AU1(x,y,t),Av(t)) I (x,y)en) defines the direction of the step. The step length e(0,1] is selected in such a way that the following function is minimized: fa =FI(U O+@A i l 0+0A FF O (u 0+C)Atiu

O

17+1A (5-10) The step direction (A~ 1

A

0 is selected in such a way that the function (5-10) becomes steadily decreasing for sufficiently small values of we(0,1]. The improved estimate is taken as the .~urent estimate of the solution of the system of nonlinear equations (5-7) corresponding to the parameter value a, and the step (5-9) I I 59 66,645 of the iterative updating scheme is repeated to determine the next improved estimate of the solution of the system of nonlinear equations corresponding to the parameter value a. This process is continued until the appropriate criteria are met, in which case the improved estimate is taken as the estimate of the velocity vector field corresponding to the parameter value a.

If the vector field (Afa,Av is not too large, then in the neighborhood of the vector field (u0,g) the nonlinear operator F° can be linearly expanded as F (u!o+AU 0fo+Afro) -P o (ro, 0,fo) +J O (il fa, Aa0, a) (5-11) Here the linear and bounded with respect to the vector field (Afi, A' 0 operator Ja o (i o ,Az, A4a) is the Jacobian of the nonlinear operator Fa(ago The operator Jag(,g 0 oAu ,A 0 a) can be more conveniently expressed in the form

J

o (a Aa Afr. -JV fr o (A u0, A f' (5-12) In the Newton method [22] the vector field (A',AV0) is defined as a solution of the system of linear equations j(U 0 (AUA, Ar) a (5-13) Then the improved estimate of the solution of the system of nonlinear equations is defined as in the relation where the scalar parameter a is taken to be equal to 1.

The improved estimate of the solution of the system of nonlinear equations which is obtained as the result of solving the system of linear equations (5-13) and then applying the relation is not necessarily a better estimate. The reason behind it comes from the fact that the Jacobian J°(fi,! 0 of the nonlinear operator Fi(i o i, 0 a is nonsymmetric and a I s~CPBIL~---BII~II~RIYI~"C ~C -L 60 66,645 ill-conditioned, so that the system of linear equations (5-13) cannot be reliably solved for the vector field The last is the consequence of the fact that the regularization constraints can only be effective in resolving the ambiguities of the system of nonlinear equations at the image points that are beyond the occluding boundaries where the quadratic form appearing in the relation (4-21) is positive definite. The regularization constraints become ineffective in resolving the ambiguities of the system of nonlinear equations (5-7) at the image points that are near the occluding boundaries where the quadratic form appearing in the relation (4-21) is not positive definite. This observation suggests that an alternative method should be employed.

In the quasi-Newton method of solving the system of nonlinear equations the system of linear equations (5-13) is replaced with the system of linear equations 4.

4 i s o a (5-14) which combines the globally convergent strategy of the well-conditioned method with the fast local strategy of the Newton method in a way that derives the benefits of both.

This quasi-Newton method defines the regularization of the system of nonlinear equations that is meaningful fox the problem of estimating the parametric velocity vector field. The linear and bounded with respect to the vector field A 0 operator (Aaa,Aia) is the approximation to the Jacobian J o o (Ai",AV) of the nonlinear operator F o defined as follows.

Let Fo(o ,goB o go be a family of nonlinear operators of the invariable copy and of the variable copy of the estimate of the velocity vector field corresponding to the parameter value a defined by the relation ~1 wv lA sraP 61 66,645 f (r 2 (p 2 +q 2 Voa 2 (r 2 (p 2+2Ve o 2 (2 2) +Yu -u 8 -fs,V s, V 0 s (a 2 +(c 2 +b 2 (S,V 2 (s,VIa) 2 (s,V 2 )1-e(a 2 +(c 2 +b (s,Vy Vla) 2) fff !fdgdxdy a, (r2+ (p 2 +q 2 JVOa 2) (g 2) 1-e pr2 +q 2 1ValoI2) (5-15) where the symbols tc, t 9,jt'o 0 are used to indicate that at every image point (x,y)en the corresponding functions depend on the invariable copy of the estimate of the velocity vector while the symbol 4t a is used to indicate that at every image point (x,y)en the corresponding function depends on the variable copy of the estimate of the Telocity vector If the variable copy a) of the estimate of the velocity vector field is identically equal to the invariable copy of the estimate of the velocity vector field, then for every BE[0,1] the nonlinear operator Fea(ia o ,ga, o ga) is identically equal to the 15 nonlinear operator Fa(ii,9°). The parameter 6 defines the degree of the feedback relaxation of the optical flow and the directional smoothness constraints through the variable copy of the estimate of the velocity vector field.

Let M(o ,vo) (AiioV be a family of.

the linear and bounded with respect to the vector field operators, where for each Oe[0,1] the operator Mea(0a,tV) (Aia,A°) is the Jacobian of the nonlinear operator Fa(Ua,a ,u r,va), considered as the function of the vector field and depending on the vector field as on the parameter, under the conditions that the vector field I 9 s tgIN I c 62 66,645 is identically equal to the vector field Then the operator Ma(fa AV 0 is defined as the member of the family '(GV, 3 V) (AdoA 0 satisfying the following requirements: 1. As a part of the globally convergent strategy, the vector field (Aii,AV") obtained as the solution of the system of linear equations (5-14) is required to be a descent direction for the functional h( o Vo,u o considered as the function of the vector field (fOgo), and depending on the vector field (iova) as on the parameter, under the conditions that the vector field be identically equal to the vector field In other words, the first-order partial with respect to the variable a derivative of the functional h( o considered as the function of the scalar variable w is required to be negative when the value of such variable o is equal to 0. It is not difficult to see that the value of this derivative is equal to _F (ca(i, a) rTo (fa, a) -lfao(o, Oa) *00 20 (5-16) The requirement for the functional (5-16) to be negative can be satisfied by in turn requiring the operator

M

0 aiO) (Ai-,A~ 0 to be positive definite.

2. As a part of the fast local strategy, the 25 operator M 0 0 (AfA',A

A

is required to be the nearest possible approximation to the Jacobian J"(V,Ge) (AWo,AVa) of the nonlinear operator FO(Oga).

The operator Mo(if 0 o 0 from the family corresponding to the parameter value 8=0.5 satisfies the 30 above requirements and is taken as the operator M(i a (Lia,AV9) appearing in the system of linear equations With this choice of the operator (AW°,A9) the system of linear equations (5-14) has the following explicit form Illlll~daaa~ 63 66,645 r a (2q aj lg o l dg G (p +q2 IV 0 2 2) 2 f(sV a (sV(Ao 0 )ds 2(a (C 2 2 b ',SV )(sV 0 2 2 Y1AQ g tOg qd g t Gr +q Vo a 1 2 (g 0 2 V VO) )dss a2+(c2+b2(s,Vr) 2 2 r 2 (gtAa dgtdxdy a fff, 2(P2 IVO 1 2 2) +y A oe r 2 (p 2 ia 2 0) 2) 2 -g -gdgtdxdy _7ov._ fff Go, r 2 2 +q 2 jVaf 2 (a 0 2 (5-17) Next, the finite-difference discretizations of the system of nonlinear equations as well as the S 5 system of linear equations (5-17) are considered. The method of constructing the estimate of the solution of the S: "system of nonlinear equations is extended to the finite-difference discretization of this system.

Although the theoretical analysis for the 10 computation of the estimate of the parametric v, .ocity vector field can be effectively performed in th, case wherr the combined generalized initial image functions rF,geG are defined over the entire set of testing functions (R 3 while the estimate of the velocity vector is defined for every image point (x,y)en and for every viewpoint teT, the actual computation can only be accomplished on finite sets of the above data sets. In this ill c-e -s Ir 64 66,645 section a finite-difference discretization of the system of nonlinear equations is described. The estimate of the parametric velocity vector field is then determined as a solution of the resulted system of a finite number of nonlinear equations.

Given a generalized function F defined on the locally convex linear topological space t(R 3 and given a testing function teP(R 3 the "convolution" of the generalized function F and the 4-esting function 0 is the infinitely differentiable func*,.u defined on the space R 3 by the relation (xy,,y ,eR (6-1) for every (x,y,t)eR.

Let F be a generalized function, and let Xel(R 3 be a given fixed non-negative testing function such go that f[X t) dxdydt=l For instance, take

R

3 x(x,y, C)0 t)/fff(9, cdyTdf(t, t,t)eR 3 where ij is

R

3 the function described later. The testing function x will S. 20 be called the "sampling function". Then the "regularization" of the generalized function F through the sampling function X is defined on the space R 3 as the infinitely differentiable function which is obtained as the convolution of the generalized function F and the sampling function X.

Let (t'k-k

K

be an increasing sequence of viewpoints containing the sequence of viewpoints k 0 Assume that for every te(t'k)k 0 the regularizations y4=(r4*X) of the generalized functions Fr assn< L.Oed with the irradiance image functions E, eS through t;i, ling 65 66,645 function X and the regularizations of the generalized functions r? associated with the feature image functions q, neH through the sampling function X are given at the points belonging to the following subset of the image plane: R h o l h e2) y) I =ilh i iEZ (6-2) Here Z is the set of integers; h'o 1 are two-dimensional real vectors; and ih'a, 1 +i 2 h'a,2 is a linear combination of such vectors with integer coefficients i l i 2 The points from the set R 2 (h' 0 2 will be called the "sampling points" corresponding to the parameter a; the functions y,Se3 will be called the "sampled irradiance image functions"; and the functions y,neH will be called the "sampled feature image functions".

The generalized functions Fr,, associated with the sampled irradiance image functions y,,ES and the generalized functions r n associated with the sampled feature image functions y,,rEH are defined on the set of 20 testing functions #(R 3 as follows. The value of the generalized function r associated with the sampled •irradiance image function y, at the testing function 0e(R 3 is given by the relation I Y y, t) x, y, t) .(xy,C)sR Q (hb x o 25 (6-3) The value of the generalized function r a associated with the sampled feature image function y, at the testing function Oe$(R 3 is given by the relation ilPil_-~ 66 66,645 t) eR 2 (hl, )x(tk) Ko (6-4) Let A (A,,XIj E,neH} be a set of real-valued constants; let g(A) be an index attached to the set 1; and let mg, {m,x,my,m ,t,m,,,mn,y,mt,,jteS,H} be a set of non-negative integer constants corresponding to the set of constants X. Then the "combined generalized sampled image function" corresponding to the index g g(A) is the 0 generalized function defined on the locally convex linear topological space D(R 3 The value of the combined generalized sampled image function rFg at the testing function oel(R 3 is given by the relation axm.xaymf,a W n ax, flay,,.Yat, r (4) 15 and is called the "observation" of the combined generalized sampled image function rF on the testing function c.

Assume the restriction to the case when each set R2(h' ,a a0, 0, an is a square grid of points on the image plane. In other words, it is a requirement that the vectors h', 1 are of equal length and are orthogonal to each other for every a Also, it is a requirement that such grids of points constitute a pyramid on the image plane which means that the following relation is satisfied:

R

2 (h h 2 )R 2 (h cR 2 (hi h (6-6) The "grid rotation method" is the sampling method that complies with the above requirements. It is commonly used to obtain numerical solutions of partial differential equations. In this method the vectors 67 66,645 h' 00 jt 1 0 2 are selected 'to be of equal length and orthogonal to each oth':-r; then -for every the vectors h Ck, 1 'hok,2 are defined as in the relations h 1, Ok ak-I 4 1 Ck.I, 2) (6-7) (6-8) It is not difficult to see that the above defined vectors haIl(T~,2are of c_,qual length and are orthogonal to each other for every and that the relation is satisfied, Let k' be a positive integer constant and *cf(R 3 be a given measurement function. Then for every off* 0 7iewpoint te~t" Ik koK and for every parameter value 15 Cra~al .an the image function g 0 and its spatial partial derivatives gc 0 (x,y,t),g~c 0 are defined on the $to":set nl as follows.

First, assume that the image point (x,y) belongs to the subset flht(hII.,,h1S 0 2 )nR 2 (h11.,i,htt, 2 Here the two-dimensional real vectors h" 11 0 2 are defined as in the :relations 0*e* (6-9) h" 0

O

2 =k1'i, 2 (6-10) while the set f1V'(h" 1 ,,h1" 0 2 )cR 2 is defined as in the relation LU~%ISldlBII3~Y)r~ *Ila~a~PI~ 68 66,645

U

N ,h y) +8,h"O,1+2h" c, 2 y) £n.

(6-11) Then for every geG the value of the component ga(x,y,t) of the image function is determined as the observation rgo(*x,y,t) of the combined generalized sampled image function F0,, on the testing function *,xy,t a defined by the relation and the values of the components of the partial derivatives 0 g"(x,y,t),gya(x,y,t) of the image function are given as the observations a xy,I) of the combined generalized sampled image function rg, o on the defined by the relations testing functions a a 1 y. Let be any image point from the set n, then for some integer numbers il,i 2 and for some real numbers 01, 02 from the interval the following relation is satisfied: y) (1z+z) (1 2 0) h",2 (6-12) Then the values of the image function g(x,y,t) and its partial derivatives ga(x,yt),gy°(x,y,t) are obtained through the bilinear interpolations g, y, (l-81) y, 1 t) +01g9 y,2a t) 2 +e (2 (X2,1y2,, t) +6,19 (x2,,y2, C) (6-13) I .g ~g eY~lirY CP q ll 69 66,645 t0 -62) (1 -01)9P (X, 1 Yl t) +6 1 gqz(X 1 2 1 t +0 2 1 1) 9- (X 2 11Y2 Vt) 0 1 gOX(X 2 2

FY

2 t) )f (6-14) g'Y(XI Y, t) (I1-0 9 g~(XI 1 y 1 3, t) +0 1 eY(XI, 21 y 1 t) (6-15) of the values of the image function go(x,y,t) and, respectively, its partial derivatives g,( 0 ,gyo(x,y,t) at the points from the set nl"(h" 0 1 1 jh".

0 2 )nR 2 0 1 ,h" 0 2 given as in the relations (XI, Iy 1 1 h"O, 1 .i 2 h 0,2' (6-16) 0 0 '2h".,21 (6-17) 0* 0* 0

(X

2 1

Y

21 I 1 1 h" 0 1 +4 (i 2 h" C 2' (6-18) 15 (X 2 2 1

Y

2 2 h" 0 (i2+ 1) 2 1 Let k" be a positive integer constant, and let h 0 ,h 0 2 be two-dimensional reO vectors defined as in the relations ho khC (6 (6-2 1) while the set fh 0 2 )c be defined as in the relation ~IBIIEl;~i~8~a~B~rPsIsll~'~ 70 66,645 0 (ha, 1 1 h 0 2

=QR

2 (ha h 0 2 (6-22) The finite-difference discretizations of the system of nonlinear equations and of the system of linear equations both corresponding to the parameter value a are defined on the set n(h, 1 ,,ho, 2 as, respectively, the system of nonlinear equations relative to the unknown estimate of the velocity vector field ,v 0 I (x,y)en(h, 1 2 corresponding to the parameter value a and its finite differences, and the system of linear equations relative to the unknown vector field I 2 corresponding to the parameter value a and its finite differences as follows.

The expressions appearing under the integral signs over the set S in the system of nonlinear equations and in the system of linear equations (5-17) have the general form Vf(x,, s) (6-23) where f(x,y,t,s) is a function of the variables (x,y)en,teT,seS, which is continuously differentiable with respect to the variables (x,y)ef for every fixed value of the variables teT,seS.

The unit circle making up the set S is replaced by a finite set, denoted by the same symbol S, of the set R 2 (ha,1,ha, 2 having the following properties: the set S does not contain the origin, and for every element s belonging to the set S the element -s also belongs to the set S. The measure ds is replaced with the point measure associating the value 0.5 to every element seS.

As a consequence of the above actions an integral over the set S of values of any expression of the form (6-23) becomes a finite sum over the set S of the values of that expression; in other words, the following relation is satisfied: -U 4 P491~-C P-AblLS~ q I I 71 66,645 f(s,Vf(x,y, ds= 0.5 (s,Vf t,s)) S seS (6-24) Given a function f(x,y,t,s) that is continuously differentiable with respect to the variables (x,y)en for every teT, seS, and given an element seS, the expression (6-23) is defined on the set R 2 in terms of the finite-difference approximations to the partial derivatives as Vf(x,y, t,s) =p,(f(x+0.5sx,y+0. 5Sy, t,s) 0 -f(x-0.5s,y-0.5sy,t,s)), (6-25) where p, is a real positive constant depending on the length of the vector s as a parameter, and the values sxsy are the components of the vector s (Sxs Taking into account the relation (6-25) and 15 the fact that the set S has the property that for every element seS the element -seS, then, in the case when f(x,y,t,s) is an odd function with respect to the variable s, the sum appearing in the relation (6-24) can be expressed as 20 pf(x+0.5sx,y+0.5syt, s).

seS S6S (6-26) It is assumed that the set G t is finite and that dg t is the point measure associating the value p(gt) to the element gteG t Here p(gt) is a real positive constant for every gteGt.

Based on the above discussion, the system of nonlinear equations can be defined on the set

R

2 2 as I P PU ~p 72 6,4 66,645 2+ (p 2 +a 2 JVC1o 2 (_a)2 Val) s~ 2

+(C

2 +b 2 (s,Vgg0) 2 (S,VgI 0 2 +Y~r(CO0-iga(0)) rE 2 2 2 Vgo 1 2 2 0 (6-27) while the system of linear equations (5-17) can be defined on the set R 2 (h, 1 ,h, 1 2 as *5 a (p2+q21VQ.12) (g-)2)2 p 3 a 2 V(A

(C

2 +b 2 (SVI_4G) 2 VG') 2 2 g~r2+(p2+q 2 IVtQI) (g 0 2 2 p V a o 0 8 W 2E 2 (S V') 2 Va.) 2 p r gt t Y, gt)XG r 2 +pq~g 2+(p 2 +q 2 VE101 2 (g")2 (6-28) The arguments of the functions appearing in the ralations 73 66,645 (6-28) have been omitted for the sake of simplicity of notation. The functions appearing under the summation over the set S sign are evaluated at the points 5sx, y+G.5sy, t) while the rest of the functions are evaluated at the points for every y) E R 2 (ho',1 1 ho, 2 teT, s=(sx, s) cS. The functions are defined, respectively, by the relations (x+A t +aI(x, Y, t) ,Y+A t'vf'(t) t+ t) gI(XA -a Y, t Lt, fc7 V t t- t+A t) (6-39) +A (xA t-GI(x,, t, y-A tV (6-30) t 2 xy t) -i 1(,y 'I(t t-z(xy t)) (6-32) The function (s,Vflo) is defined by the relations Yq~L~ls~L~aPumuna~i~s~l~sa~llsl---~ 74 66,645 Vgu)) Vu° 5s,y+0. 5sy, t) P(s u(x++Sx,y+sy, t) (6-33) while the function is defined as in the relations (s,V(A oo)) -(s,V(AS (x+0.5 9x,y+0.5sy, P,(Au (x+Sx, Y+Sy, t) -AQu(x,y, (6-34) The function b 2 (s,V'gt0) 2 is defined by the relation b2 V 0r) 2 p(g) b V 2 (6-35) where each function (s,V'Ia) 2 ,gtG t is determined as 0(sV~~ (S,V 5sy+0. 5sy, t, o 2 10 V'g? 2+ Vg y+sy, t, CO, 2.

(6-36) Above, for each value (equal either to or to (x+sx,y+sy,t)) the value of the function (s,V'gtO(x',y',t,1Vg))2 is given by the relation The restriction of the system of nonlinear equations (6-27) and the restriction of the system of linear equations (6-28) on the set n(h,,l,ha, 2 of the set R 2 1 2 are both accomplished through introduction of the boundary condition as follows. For every image point (x,y)en(ho,l,ho, 2 and for every element s (Sx,Sy)ES such that the point is outside of the set n(h,,,ho, 2 the value of the estimate of the velocity vector (G°(x+SxY+ y,t) is defined to be identically equal to the value of the estimate of the velocity vector Under these conditions the system of nonlinear equations (6-27) and the system of linear equations (6-28) are well defined for every image point from the set n(h, 1 ,h, 2 and, therefore, form, i, r-U E ~CB ~PI~P- I C 75 66,645 respectively, the system of nonlinear equations relative to the unknown estimate of the velocity vector field (ua(x,y,t) I (x,y)en(h, 1 ,h o 2 corresponding to the parameter value a and its finite differences, and the system of linear equations relative to the unknown vector field (Au(x,y,t) ,Ag I corresponding to the parameter value a and its finite differences.

Let us denote by the symbol ISI the number of elements of the set S, then for every s (sx,sy)eS the constant Ps can be defined as in the relation p=(s (6-37) The following are two commonly used ways of selecting the set S S-hoa, 1 -he,l, h, 2 2 (6-38) Saha,,, -h 0 11 h 0 0 2 ,h 0 1 0 2' s-a, o 2' -ho, 2' h, I+ha 2 -he, 1 -ha, 2 a 2 -ha,i +ha, 2f.

(6-39) If the set S is selected as in the relation (6-38) we say that the "five-points finite-difference discretization" is used, while for the selection specified by the relation (6-39) we say that the "nine-points finite-difference discretization" is used.

Let k and let y, t) (h k, hok,)) be an estimate of the velocity vector field corresponding to the parameter value Uk. Then the initial estimate y, t) V'1 y) cQ (hak,, of the velocity vector field corresponding to the parameter value ak+I is kP~I b~4~P~~ 76 66, defined as follows. For every y) eCQ (hak,,h Ck1, such that eQ (h a 0 1 1j, 2 the value of the estimate a ak-t(t.) is given by the relation while for every (x,y)eQ(h Ok~l hak12 such that y) 0(hal,4 hak 2 the value of the estimate (a Ok1 t) Vk-t( t) is given as in the relation 645 a 0*

VS

*1 S a

S.

0 V. V C S *5 a

V

.VVG

V

It, Ok (X V kf 10 (6-41) where (6-42" The system of linear equations (6-28) is symmetric and positive definite. As shown below, the sy ,tte of linear equations (6-28) can be partitioned into two systems, one of which can be directly inverted. Based on these properties a fast preconditioned iterative niathod is.

constructed. If vwe introduce an orderi4ng on the set r1(h, 11 h, 12 then the system of linear equations (6-28) can be expressed in the matrix notation as ze 00, VO) 6L ill, AV*) 1= -r VO) 77 66,645 where (AQ, OF)' is the transpose vector of the row-vector (Ai, Ay') corresponding to the vector field ({(Afig y,t) ACP(t) I y) n(h,, 1 2 I which is obtained by first listing the elements of the set A y, t) I y) s n(h, 1 2 according to the specified ordering, then listing the element A1~ 0 The "natural ordering" on the set 2 can be defined as follows.

The relations (6 and (6-22) imply that every element (x,y)efl(h,, 1 ,h 0 2 can be uniquely represented in the form _=ihai,+' 2 ha 2 where and i 2 are integer numbers.

Then the element ai 1 h, 0 1 +i 2 h, 0 2 6f(h, 3 1 1 2 follows the Aft. element Ei i' 1 ho 0 1 2 h,, 2 Ef(h,, if i 2 >i' 2 or if i 2 i' 2 and i 1 >i' 1 Let the matrix be the diagonal part of the matrix while the matrix 0 )is the off-diagonal part of the matrix M 0 i 0 ~,so that the ~.~,.relation ft (7-2) 20 is satisfied, then the system of linear equations (6-28) takes the form ft (7-3) 0 6 Let C 0

(V

0 be a diagonal matrix defined by 25 the relations V) 0 V)12 (7-4) then the following relation is satisfied Let (Au',Avl) be a vector defined by the relation (7-6) 78 66,645 If both sides of the system of linear equations of the unknown vector (Afi,A°) are multiplied on the matrix Co(~goa)T from the left and if the unknown vector (AW",AVi) is stituted with the unknown vector the system of linear equations becomes the following system of linear equations relative to the unknown vector (AQuA') A T=Cc VO) -TB1 r) c 0 (a o V) 1

A

b T -c a e) -Tf (s a t) (7-7) Let be a matrix defined as in the relation V) V)

T

BO(, V

O

-r (7-8) while h°(ia°,V 0 is a vector given by the relation

O(

O

V) =cr VO) -TBp ,fi 15 (7-9) then the following basic iterative method can be used to obtain an approximation to the solution of the system of linear equations (6-28).

*'The initial approximation (AIo',AVo) to the solution of the system of linear equations is defined to be identically equal to zero. For every the approximation (AUn+,,On1 0 is defined in terms of the approximation (AanO,AVn 0 as in the relation o (A il^ G (AiiuAiL)

T

-ha (i (7-10) The process is continued until a proper approximation (AUN,, AYa to the solution of the system of linear equations is achieved. The approximation (Aio,A) to the solution of the system of linear equations (6-28) is defined by the relation 79 56,645 (A A

T

=CO jV) -1 7AN) T~ (7-11) The performance of the basic iterative method (7-10) can be improved with the help of the polynomial acceleration applied to the basic iterative method. For the sake of simplicity of the presentation the following notation shall be used: w= (A il, As) G= G I, f) ,hbeh llf (7-12) Then the basic iterative method (7-10) becomes w. =Gw,-h.

(7-13) s. The conjugate gradient polynomial acceleration of the basic iterative method (7-13) can be described as 15 follows. Start with the initial approximation bold w 0 which is identically equal to zero, then for every the following iterative procedure is applied:

S

1 ="Pn (Yn (Gn-h) (l-Yn) W n (l-Pn) W-3.

(7-14) i 20 Here the coefficients p, and n, are given by the relation n- 1+ ann nn n n-1 Pn 1po=1+ Y an-an n (7-15) where an- p, ,n 0, n= 0=0, n 7) nl,-=O (7-16) pn=r+Pnpn-_Ln>lPo=ro, q=Pn -QGp,naO, (7-17) s -uesl--r I I~OP~ D ls 80 66,645 r,=Gvw-w-hb,n>O.

(7-18) The conjugate gradient method can be expressed in the alternative form 1 gg (7-19) where the coefficients an are defined as in the relation while the vectors pn, qn are given by the relation (7-17).

In an exemplary embodiment of the invention, 0 the five level resolution pyramid was used with the value of the parameter a decreasing by a factor of 2.0 for each successively finer resolution level. A nine-points finite-difference discretization was used. The positive 15 integer constant K' appearing in the relations (6-10) is equal to 4, while the positive integer constant K" appearing in the relations (6-21) is equal to 2.

With these values, the present invention is capable of estimating the velocity vector field with good accuracy in the case of the images of real world scenes containing multiple occluding boundaries.

The many features and advantages of the P invention are apparent from the detailed specification and thus it is intended by the appended claims to cover all such features and advantages of the invention which fall within the true spirit and scope of the invention. Further, since numerous modifications and changes will readily occur to those skilled in the art, it is not desired to limit the invention to the exact construction and operation illustrated and described, and accordingly all suitable modifications and equivalents may be resorted to, falling within the scope of the invention.

-81- Parts List: Camera(s) 12 Digitizer 14 Computer 16 Display device Pixel 22 Intermediate image 24 Image 26 Image 28 Vector or line Two-dimensional area 32 Two-dimensional area 38 Horizontal axis 50 Right digital image 52 Left digital image 54 Establish correspondence between left and right images 56 Estimate intermediate disparity vector field 58 Estimate intermediate disparity vector field Interpolate and interlace of intermediate images 62 Depth image Select initial estimates of intermediate disparity vector fields 72 Compute values of the filtered image band and partial derivative of the filtered image band 74 Take initial estimates of intermediate disparity vector fields 76 Form the system of linear equations 78 Solve the system of linear equation III l -~Pe 1LT- I -82- Parts List Continued: Current estimate of the intermediate disparity vector field improved to a desired degree 82 Current level of resolution finest 84 Increment the current level of the multilevel resolution pyramid to the next finer level Project the resultant estimate of the intermediate disparity vector field to the initial estimate of the current finer level of resolution e *oo• 9 o oO a o *oo I II Ir I-'-411(IPC~ I_ _I -83- Apendix A--Docket 66,645 S1993 Eastman Kodak Company, Rochester, N.Y. 14650-2201 METHOD AND APPARATUS FOR CONSTRUCTING INTERMEDIATE IMAGES FOR A DEPTH IMAGE FROM STEREO IMAGES

BY:

Sergei Fogel ac a.

a 0 1993 Eastman Kodak Company A portion of the disclosure of this patent document contains material which is subject to copyright protection.

The copyright owner has no objection to the facsimile reproduction by anyone of the patent document or the patent disclosure as it appears in the'Patent and Trademark Office patent file or records but otherwise reserves all copyright rights whatsoever.

a Staas Halsey 1825 K Street, N.W.

Washington, D.C. 20006 202-872-0123 i -rpasl- r sR -84- #define OSIZE #define ASIZE #define BANDS 3 #define BUTFEL 4000 #define SMTH 2.0 2.0 #define STPD #define STEP 0.4 #define SGMB 2.25 #define SGMQ 2.25 #define SGMW 2.25 #define RHG 1.0 0.5 #define R.HX 3.0 0.5 #define DSG 1.0 1.0 #define MSG 2.0 1.414213562 #include .cmath~h> #include <stdio.h> #include <sys/file~h> #include "iopackage.h" *#include "cor _corslv.h' struct opt-flow (long kuvrn[ASIZE], kuvp[ASIZE]; float *ggm[BANDS], *gxm[BMNDS], *gym~[BANDS], *sxm[BANDSI, dhm[ASIZE]; float *xgm[BANDS], *xxm[BANDS], *xym[BANDS], *sym[BANDS], dvmIIASIZE]; float *ggp[BANDS], *gxp[BM.NDS], *gyp[BANDS], *sxp[BANDS], dhp[ASIZE]; float *xgp[BANDS], *xxp[BANT)S], *xyp[BANTJS], *syp[BANDS], dvp[ASIZEI; float *fgm, *fgp dtm[(ASIZE], dtc[ASIZEI, dtp[ASIZE];1 optf; struct ars..data (long cns, cnv[ASIZE]; long *mll[ASIZE], *m12[ASIZE],*m~1[ASZEI, *w.22[ASIZE]; float *wll[ASIZE], *wl2[ASI, *w21[ASJZEI, *w22[ASIZE]; float wmn[ASIZE], wxyIIASIZE], wuv[ASIZE], dlt[ASIZEI; float *auul, *auu2, *auu3, *auu4, *auu5, *ainv, *bvvl; Sfloat *cuvl, *g2sq, *g3sq, *g4sq, *g5sq, *dusq, *bjnv; float *pjcvi, *pjjcv, *qucv, *zjyjc., a, b, bsq, c, csq; float *rcv *pv7v, *qvcv, *zvcv, r, p, psq, q, qsq; __float *uuvf, *u~vf, *udvf, *utvf, *duvf; Wifloat *vvvf, *v(Jvf, *vdvf, *vtvf, *dvv,r arsd[ASIZE]; double sgms, sgds, sgm 4.0, sgmb, scln, tmsg MSG, scl; double sgmg, sgdg, sgd 2.0, sgmq, sgmw, dsg DSG, two float *hv~h0, *wrk00, *wrk0l, *wrk02, *wrk03, *wrk04, float *vvOhO, *wrk06, *wri~k07, *wrqk(J8, *wrk09, *wkO *wrkllI; float *wv .kl2, *wrk3, *wrk14, *wrk15, *wrik16; long *mv~dhO, *mvphp, *mv~hp, *mvmJhp, *mvrnh0, *mvrrim; long *mv~hm, *mvphn, *mvph0, *mllp, *mp22p; long *mllm, *m]2rn, *m.p2lm, *mp22m, *ml2p, *w2lp; float gun 0.00 1, gul 0.0 1, rhg RIIG, isc 0.05; float gvn 0.001, gvl 3.01, rhx RHX, thre; long k-vOhO, kvphp, kvOhp, kvmhp, kvmh0, kvrnhm, kv~hm, kvphm, kvphO; long di 0, rho 0, msp 0, lvs 0, 1ev, ksp, nij; long dkj 0, tau 0, nsp 0, 1st 0, nhv, nxy, bns 0; /~The program. for estimating the velocity vector fields /~from a given time-varying image sequence. main(argc,argv) 1* 1* nsp the variable of the integer type, which is equal to 1* the number of vector fields that are involved in the 1* computation. 1* vs the variable of the integer type, which is equal to 1* the number of resolution levels that are used in the computation. st the variable of the integer type, which is equal to 1* the resolution level at which the estimation process 1* stops. /~rho the variable of the integer type, which is equal to 1* the degree of subsampling of the velosity vector fields on the finest level of the resolution piranid relative to subsaxnpling of the image functions and their partial derivatives. /~tau the variable of the integer type, which is equal to 1* the degree of subsampling of the image functions and 1* their partial derivatives relative to the sampling of 1* the initial image functions.

/~ust the variable of the float type, which is equal to the 1* x-component of the initial estimate that is used in the 1* computation. 1* vst the variable of the float type, which is equal to the 1* y-component of the initial estimate that is used in the computation. nh the variable of the integer type, which is equal to the number of vectors along the horizontal direction that are involved in the computation. nv the variable of the integer type, which is equal to /*the number of vectors along the vertical direction that are involved in the computation. I* kx thevarible of the integer tp, which is equal to the number of pixels defining the left and the right /*offsets for the image functions and their partial 1* derivatives; the number kx must be a multiple of the number two in the power of. lvs in the case when the number lvs is odd, and lvs I in the case when the number lvs is even. ky the variable of the integer type, which is equal to the number of pixels defining the lower and the upper offsets for the image functions and their partial 1* derivatives; the number ky must be a multiple of the number two in the power of: lvs in the case when the number lvs is odd, and lvs I in the case when the number lvs is even. nx the variable of the integer type, which is equal to the number of vectors along the horizontal direction. /~ny the variable of the integer type, which is equal to 1* the number of vectors along the vertical direction. /~sgmn the variable of the type double specifying the value of the horizontal component of the sigmaa factor. 1* sgd the variable of the type double specifying the value of the vertical component of the sigma factor.

/~msg the variable of the type double specifying the b~Be~lP~111188PI ~C -86multiple change factor for the value of the sigma factor sgm between the resolutions. dsg the variable of the type double specifying the multiple change factor for the value of the sigma S* factor sgd between the resolutions. isc the variable of the float type specifying the initial image scale factor. rhg the variable of the float type specifying the weight for the optical flow constraints corresponding to the image functions.

rhx the variable of the float type specifying the weight for the optical flow constraints corresponding to the x-directional partial derivatives of image functions. gun the variable of the float type specifying the weight for the u-regularization constraints appearing in the system of nonlinear equations. gvn the variable of the float type specifying the weight for the v-regularization constraints appearing in the system of nonlinear equations. I* gul the variable of the float type specifying the weight for the u-regularization constraints appearing in the system of linear equations. gvl the variable of the float type specifying the weight P for the v-regularization constraints appearing in the system of linear equations. rep the variable of the float type defining the criterion S for termination of nonlinear iterations. fep the variable of the float type defining the criterion for termination of linear iterations. mxn the variable of the integer type specifying a maximal number of nonlinear iterations. mxu the variable of the integer type specifying a maximal number of trials allowed in the attempts of updating a nonlinear iteration. S mxl the variable of the integer type specifying a maximal number of linear iterations. mxs the variable of the integer type specifying a maximal number of times the linear system solver can be used on each nonlinear iteration. intargc; char *argv[]; unsigned char; char *fldt, *flmc, *flpc; double sqrt0, dtmp, stmp; long i, n, nx, ki, ni, dd, mh, itst, nh 0, kx 0; long j, k, ny, kj, nj, dn, mv, intv, nv 0, ky 0; long mxn 0, mxl 0, imc 0, jmc 0, fmc 0, fdt 0; long mxu 0, mxs 0, ipc 0, jpc 0, fpc 0, ckk, bn; float avrgm[BANDS], *ggg, *ggx, *ggy, *gsx, tm, ust 0.0, rep 0.00001, *fgg, tp; float avrgp[BANDS], *gxg, *gxx, *gxy, *gsy, tO, vst 0.0, fep 0.00001, temp; for (k 1; k argc; if -87if Vd) (*(argv[k] 44) 'T) fldt argv[k]+5; fdt+-i; if fdt 1) fprintf(stderr,'Too many fldt files.\n"); exit(1); }else if 'in) (*(argv[kI+4) V) flmc argv[k]+5; fmnc++; if fmc 1) fprintf(stderr,'Too many flinc filesM' t exit( 1); AOL else if V) flpc argv[k]+5; fpc++; if (fpc 1 fprintf(stderr,"Too many flpc files.\n"); exit(l); else fprintf(stderr, 'W)7The parameter %s" t argv[k]); fprintf(stderr, is not recognizable.\n"); exit(l); else if if (*(argvfk]+2) 'in) V if (sscanf(argv[kJ, &imc) 1) fprintf(stderr,'WO7Can't read argv[k]); fprintf(stderr," it must be -inc,<long>\n"); exit(l); else if (*(argvk]+2) V if sscanf(argv[k], &ipc) 1) fprintf(stderr,"\O7Can't read argv[k]); fprintf(stde-r," it must be -ipc~long>\n"); exit(l); Ielse if T s) (*(argvfk]+3) =V -88fprintf(stderr," it must be -nv<long>\n"); exit(1); else if TI) if (*(argvMk+3) =T if sscanf(argvllk], &lvs) 1) fprintf(stdlerr,'"OOCan't read argv[k]); fprintf(stderr," it must be -lvs<long>\n"); exit(l); Ielse if Ts) (*(argvk]+3) =T if sscanf(argv[k], &lst) 1) fprintf(stderr,"\O7Caxi't read argv[k]); fprintf(stderr," it must be -Ist'dong>\,n"); exit(l); if (*(rgvkJ)= ==T if (*(ajrgv[k]+2) (*(airgv[k]+3) V if sscanf(argv[k], &tau) 1) fprintf(stderr,"\OO7Can't read argv[k]); fprintf(stderr," it must be -tau<long>\n tt exit(1); else if (*(argy[kI+1) V if (*(ajrgv[k]+2) =T if sscanf(argv[k], &ust) 1) fprintf(stderr,'VO7Can't read argv[kI); fprintf(stderr," it must be -utfoa>n) exit( 1); Ielse if V) if (*(argv4kI+2) Ts) (*(airgv[k]+3) =T if sscanf(argv[k], "-vst%f, &vst) 1) -89fprintf(stderr,'4)7Can'It read argv[k]); fprintf(stderr," it must be -vst~float>\n"); exit(l); Ielse if (*(argvk]) (*(argv[kI+1) Yk)) if *(argvMk+2) Y x' if sscanf(argv[k], "-kx%d t &kx) 1) fprintf(stderr,'NDO7Can't read argv[k]); fprirnf(stden-," it must be -kx<long>\n"); exit(l); Ielse if *(argv[k]+2) if sscanf(argv[k], &ky) 1) fprintf(stderr,N)O7Can't read ar~v[k]); *fee* fprintf(stderr," it must be -ky<long>\n~"); Poe*exit( 1); Ielse if if (*(airgv[k]+2) Wh) if sscanf(argv[k], &rho) 1) fprintf(stderr,'\OO7Can't read argv[k]); fprintf(stderr," it must be -rho~long>\nI); exit(1); .6.6 )else if (*(argv[kI+2) Wh) if sscanf(argv[k], &rhg) 1) fprintf(stderr,'\007Gan't read .argv~k]); fprintf(stderr," it must be -rhg<float>\n"); exit(1); Ielse if Wh) (*(argv[kI+3) Y if sscanf(argv[k], t -rhx%f', &rhx) 1) fprintf(stderr,'V\O7Can't read argv[k]); fprintf(stderr," it must be -rhx<float>\n'); exit(l); )else if (*(argv~kI) (*(argvk]+l) if Vu) W if sscanf(argv~k], &guni)! I1 fprintf(stderr,N)O7Can't red T ,,arv fprintf(stderr," it must be -gun<rioacAn"); exit(l); Ielse if if sscanf(argv[k], &gul) 1) fprintf(stderr,'NW7Can't read argv[k]); fprintf(stderr," it must be -gukfloat>'n"); exit( 1); Ielse *if (*(argv[kI+3) **if sscanf(argv[k], "-gvn%f &gvn) I 1 fprintf(stderr,'VO7Can't read argv[k]); fprintf(stderr," it must be -gvn<float>\n"); exit( 1); Ielse if (*(airgv~kI+2) T1))

I

if sscanf(argv[k], &gvl) 1) 0 fprintf(stderr,'N)7Can't read argv[k]); fprintf(stderr," it must be -gvl<float>\n"); exit(1); Ielse if W if (*(argv~k1+2) Yx) V if sscanf(arkv[k], &mxn) I) fprintf(stderr,'NO7Can't read argv[k]); fprintf(stderr," it must be -nixn<long>\n"); exit(1); I else "*(argv[kI+2) Yx) V -91if sscanf(argv[k], "-mxu%1i", &rnxu)!=1 fprintf(stderr,'WOCan't read argv[k]); fprintf(stderr," it must be -mxu'clong>\n"); exit(1); else if (*(arg[klJ+2) Yx) T'1)) if sscanf(argv[k], "-mxl%d' t &mxl) 1) fprinrf(stderr,N)O7Can't read argv[k]); fprintf(stder-r," it must be -mxklong>\n"); exit( 1); else if Yx) T) if sscanf(argv[k], '"4nxsOd", &mxs) 1= I) fprintf(stderr,'V)O7Can't read %sIt, argv[k]); fprintf(stderr," it must be -mxs<long>\,n"); exit(l); if I fpit*sder "fd **nfd) pitfsdu,"itnfd) if ((fmc= fprintf(stderr,"h "flm e %,mc isno %dpjc lcic;m) fprintf(stderr,"e "fl e is no %dspeci d\n'")lpic;p) fprintf(stdut, t flparamets r ipc eihe not bee&n, sep, pcj); fprintf(stderr,"Trhasfilen ms o thwrn spe lufe ns "-92while scanf("%d', &nsp) 1) fprintf(stderr,"\O7Can't read read input for fprintf(stderr,"nsp please try again; nsp if( nsp >ASIZE) fprintf(stderr,'The parameter nsp, has been set to the value\n"),; fprintf(stderr, "that is greater than the upper limit of the~n"); fprintf(stderr, "size ASIZE of array of structures arsd; nsp= while scanf("%d", &nsp) 1) fprintf(stderr,N)O7Can't read read input for') fprintf(stderr,"nsp please try again; nsp while (lvs <1) fprintf(stderr,"The parameter lvs either has not been seM"); fprintf(stderr,"or has been set to the wrong number lvs= while scanf("%d", &lvs) 1) fPrintf(stderr,'WO7Can't read read input for") fprintf(stderr,"lvs please try again; Ivs while rho 1) fprintf(stderr,"The parameter rho either has not been set\n"); printf(stderr,"or h.qs~ been set to the wrong number; rho= while scanf("%d", &rho) 1) fprintf(stderr,"\OO7Can't read read input for") fprintf(stderr,"rho please try again; rho ::while (tau <1I fprintf(stderr,"The parameter tau either has not been setWn); fprintf(stderr,'or has been set to the wrong number, tau while scanf("%d", &Mau) 1) fprintf(stderr,'WO7Can't read read input for fprintf(stderr,'tau please try again; tau while nh 1) fprintf(stderr,"The parameter nh either has not been set\n"); fprintf(stderr,"or has been set to the wrong number; nh= while scanif(%d", &nh) 1) fprinf(stderr, M.O7Can't read read in~put for -93fprintf(stdcrr,"nh please try again; nh while (fly 1) fprintf(stderr,"The parameter nv either has not been set\n'); fprintf(stderr,"or has been set to the wrong number-, nv while scanf( 1 &nv) 1) fprintf(stderr,'N0O7Can't read read input for fprintf(stderr,' t nv please try again; nv' for k 1; k <argc; if if (*(argvfk]+2) (*(argvk]+3) if sscanf(argv~k], &ksp) 1= I) fprintf(stderr, "\O7Can't read argvjjk]); fprintf~sderr," it must be -ksp<long>\,n"); exit(l); )else if((ksp<O0) 11(ksp >nsp) fprintf(stderr,"MO7ksP %d is set to wrong number.\n, ksp); exit(l); jf (*(argv[k+l) f((*(airgv[k]+2) W) S. if sscanf(argv[k], &intv) 1= 1) fprintF(stderr,N)O07Can't read argv(k]); fprintf(stderr,", it must be: -cns<long>\n 9 exit( 1); Ielse if (intv 0) (intv nsp (ksp nsp) arsdfksp].cns =intv; Ielse if (ksp nsp) fprintf(stderr,"OO7ksp %d is set to wrong number.\n", ksp); exit(l); Ielse if ((intv 0) 11 intv >=nsp)

I

fprintf~stderr,"\OO7Can't read cnis", argv[k]); fprintf(stderr,', is set to the. wrong numnber.An"); -94exit( 1); for( k= 1; k<argc; if T if (*(ajrgv[k]+2) T s) (*(argv[kI+3) if sscanf(argv[k], &ksp) =1 fprintf(stderr,"'OO7Can't read argv[k]); fprintf(stderr," it must be -ksp<long>\n"); exit(l); else if ksp 0 )11 ksp nsp) 0 fprintf(stderr,'V\07ksp %d is set to wrong number.\n", ksp); exit(l); ifls if 0 if sscanf(argv[k], &temp) 1) fpriinf(stderr,NO7Can't read %s it argv[k]); fprintf(stderr, "must be: -asp<float>\n");, exit(1); 0.0 if((ksp (ksp <nsp)) 0 arsdfkspl.a temp; Ielse if (ksp nsp) for ksp 0; ksp nsp: ksp++)

I

a 3dfksp].a trmp; ksp nsp; Ielse if (*,agv4k]) if Ts) (*(argv[kI+3) =p' if sscanf(argv[k], &ternp) 1) fprintf(stderr,"\O7Can't read %s it argv[kJ); fprintf(stderr,"must be: -csp<float>Nn'); exit(l); if (ksp 0) (ksp <nsp)) arsdfkspl.c temp; else if ksp nsp) for ksp 0; ksp nsp; ksp-i+) arsdfkspl.c temnp; ksp nsp; Ielse if Yk) V) if sscanf(argv[k], 1) fprintf(stderr,W~7Can't read argv[k]); fprintf(stdlerr,", it must be -ckk<long>\n"); exit( 1); Ielse if (ksp nsp) fprintf(stderr,'\0O7ksp %d is set to wrong number.\n", ksp); exit(1); )else if((ckk <0 11(ckk> arsd[ksp].cns) fprintf(stderr,'\007ckk %d is set to wronig number.\", ckk); exit(l); Ielse if (*(rgvfk+2) (*fargv[k]+3) if sscanf(argv[k], &intv) 1) a:..:fprintf(stderr,'%0O7Can,'c revid 9/s it must be", argv[ki); fprintf(stderr,": -cnv<Itong~n"); exit(l); else if ksp nsp) fprintf(stderr,N)O7ksp %d is set to wrong numberNn", ksp); exit( 1); arsdilkspl.cnv[ckk] intv; else if (*(argv4~) (*(aIrgv[k]+1) V if (*(argvkI+3) =p'

I

-96fprintf(stderr,'WO7Can't read %s it argv[k]); fprintf(stderr,"must be: -bsp~floarAn' t exit(1);' if ksp 0 ksp nsp) arsd[ksp].b temp; else if ksp nsp) for ksp 0; ksp nsp; ksp++) arsd[ksp].b temp; ksp nsp; else if (argv[k]) (agv Y) Ts) 'p' if sscanf(argv[k], &temp) 1) fprintf(stderr,NX)7Can't read %s it argv[k]); too: fprintf(stderr, "must be: -rsp~float>\n"); if(ksp 0 ksp nsp) arsd[ksp].r temp; Iels.-if (ksp==nsp) :for (ksp=O0; ksp< asp; ksp++) arsd[ksp].r temp; Wit ksp =nsp; else if 'p' if (*(argv[kI+2) Ts) 'p' if sscanf(argv[k], &temp) 1) fprintf(stderr,"'O7Can't read %s it argv[k]); fprintf(stderr,"must be: -psp<float Wn); exit(l); if C((ksp 0) (ksp nsp)) arsd[ksp].p teinp; Ielse if ksp nsp) -97for ksp 0; ksp nsp; ksp-Harsd[ksp].p temnp; ksp nsp; else if (*(airgv[kI+l) if Ts) =p' if sscanf(argv[k], &temp) 1) fprintf(stderr,"\V07Can't read %s it argv[k]); fprintf(stderr, "must be: -qsp~floatz>n"); exit(l); Aft if (ksp 0) (ksp nsp)) arsd[ksplj.q temp; )else if (ksp=--nsp) for (ksp ksp <nsp; ksp++) arsd[ksp).q =temp; ksp sp; I else if W) if W h) ~if sscanf(argv[k], "-dhm%f &temp) =1 fprintf(stderr,NL)7Can't read %s it must be", argv[k]); fprintf(stderr,": -dhm<float>\n"); exit( 1); if C (ksp 0) (ksp nsp)) optf.dhmnhksp] temnp; Ielse if ksp nsp) for (ksp=0; ksp <nsp; ksp++) optf.dhm[ksp] temp; ksp nsp; Ielse if Wh) 'p' -98if sscanf(argv[k], &temp) 1) fprintf(stderr,"\O7Can't read %s it must be", argv[k]); fprintf(stderr,": -dhp<float>\n"); exit(1); if ksp 0) (ksp nsp)) optf.dhp[ksp] teinp; else if ksp nsp) for ksp 0; ksp nsp; ksp+-ioptf.dhp[ksp] temnp; ksp nsp; else if 'mn)) *if sscanf(argv~k], &temp) 1) **.*fprintf(stderr,"07Can't read %s it must be", argv~k]); 0. fprintf(stderr,": -dvm.~float>\n"); exit(l); if ((ksp 0 ksp nsp) optf.dvin[ksp] temp; else if (ksp nsp) :for ksp ksp nsp; ksp++) optf.dvm[ksp] temp; Iksp =nsp; Ielse if (*(argv[kI+2) V) if sscanf(argv~k], "-dvp%f &temp) =1) fprintf(stderr,"'W7Can't read %s it must be", argv[k]); fprintf(stderr,": -dvp<float>\n"); exit( 1); if (ksp 0) (ksp nsp)) optf.dvp[ksp] temp; Ielse if ksp nsp) for ks-p 0; ksp nsp; ksp++)

I

optf.dvp[ksp] temp; -99ksp nsp; Ielse if (*(argv[kI+3) if sscanf(argvfk], &ternp) 1) fprintf(stderr,"\O7Can't read %s it must be", argv[k]); exit( 1); if ksp ksp nsp) optf.dtm[ks-p] =temp; Ielse if ksp nsp) for (ksp ksp nsp; ksp++) optf.dtm[ksp] =temp; ksp =nsp; Ielse if (*(ajrgv[k]+3) V) if sscanf(argv[k], &temp) 1) fprintf(stderr,'WO7Can't read %s it must be", argv[kD); fprintf(stderrn": -dtc<float>\n"); exit(l); if (kp (ksp <nsp)) optf.dtc~ksp] temp; Ielse if (ksp nsp) *for (ksp ksp <nsp; kspi+) optf.dtc[ksp] temp; ksp =nsp; Ielse if 7 if sscanf(argv[k], &temp) 1) fprintf(stderr,NO7 Can't read %s it must be", argv[kJ); tbprintf(stdex.": -dtp<float>Nn"); exit(l); if (ksp 0) (ksp nsp)) optf.dtp[ksp] temp; Ielse if ksp nsp -100for ksp 0; ksp nsp; ksp++) optf.dtp[ksp] temp; ksp nsp; Jelse if T) if sscanf(argv[k], &temnp) 1) fprintf(stderr,'VX7Can't read %s it must be", argv[k]); fprintf(stderr,": -dlt<float>\n"); exit(l); if (ksp 0) (ksp nsp)) arsd[ksp].dlt[ckk] temp; )else if ksp nsp) for (ksp ksp <nsp; ksp++) arsd[ksp].dlt[ckk] temp; a 0 ksp =nsp; else if (argv[k]) =w **if

W)

if sscanf(argv[k], "-winn%f', &temp) 1) fprintf(stderr,N)O7Can't read %s it must be", argv[k]); fprintf(stderr,": -wmn<float>\n"); exit( 1); if (ksp (ksp nsp)) arsdllksp].wmn[ckk] temp; else if ksp nsp) for (ksp ksp nsp; ksp++) arsd[ksp].wmn[ckk] temp; ksp risp; Ielse if Yx) 'y' if sscanf(argv[k], &temp) 1) -101fprintf(stderr,"'VO7Can't read %s it must be", argv[k.]); fprintf(stderr,": -wxy<float>\n"); exit(1); if (ksp 0) &&(ksp nsp)) arsd[ksp].wxyfckkl temp; Ielse if ksp nsp) f for ksp 0; ksp nsp; ksp++)

I

arsd[ksp].wxy[ckk] temnp; ksp =nsp; Ielse if (*(airgv[k]+2) Vu) if sscanf(argvlk], "-wuv%f &temp) =1 fprintf(stderr,'W)0Can't read %s it must be", argv[k]); fprintf(stderr,": -wuv<float>\n"); exit( 1); *0if ((ksp &&(ksp <nsp)) arsd[ksp].wuv[ckk] =temnp; Ielse if (ksp nsp) for ksp 0; ksp nsp; ksp++)

I

arsd[ksp].wuv[ckk] temnp; ksp nsp; fprinrf(stderr, "lys rho lvs, rho); fprintf(stdout, "lvs rho lvs, rho); fprintf(stderr, "nh kx rnxn mxu nh, kx, mxn, mxu); fprintf(stdout, "nh kx mxn mxu dn", nh, kx, mxn, mxu); fprintf(stderr, "nsp tau ",nsp, tau); fpi-intf(stdout, "nsp tau ",nsp, tau); fprintf(stderr, "nv ky mxI mxs W dn", nv, ky, mxl, mxs); fprintf(stdout, "nv ky mxl mxs dn", nv, Icy, mxl, mxs); fprintf(stderr, "ust %8.5f, rhg ",ust, rhg); fprintf(stdout, 'lust rhg ",ust, rhg); fprintf(stderr, 'rhx gun %8.5f, gui. rhx, gun, gui); fprintf(stdout, "rhx gun %8.5f, gui rhx, gun, gui); fprintf(stderr, "vst %8.5f, isc %8.5f, ",vst, isc); fprintf(stdout, "vst isc ",vst, isc)!; fprintf(stderr, "gvn %8.5f, gvI= %8.56an", gvn, p1l); fprintf(stdout, "gvn %8.5f, gvI gvn, gvi); for ksp 0; ksp nsp; ksp++ -102- (pit~teT "s d s) fprintf(stderr, "ksp ksp); fprintf(srdout, "as arksp a fprintf(stderr, "a %6.3f; ",arsd[ksp].a); fprintf(srdout, "ac %6.3f; arsd[ksp.a); fprintf(stderr, "c arsd[ksp].c); fprintf(stdout, "b %6.3f; arsd[ksp].c); fprintf(stderr, "b arsd[kspl.b); fprintf(stdout, "r %6.3f; arsd[ksp].b); fprintf(stderr, "r arsd[ksp].r); fprintf(stdeut, "rp %6.3f; arsd[ksp].r); fprintf(stderr, "p arsd[ksp].p); fprintf(stdeut, "p =r %6.3f; arsd~ksp].p); fprintf(stderr. "q arsd[ksp].q); fprintf(stdout, "q arsd[ksp].q); fprintf(stderr, "cns arsd[ksp].cns); for ksp 0; ksp nsp; ksp+-i- (pit~ter "s d s) fprintf(stderr, "ksp ",ksp); fprintf(stdout, "ksp kspdh~sp fprintf(stderr, "dhm optf.dhm[ksp]); fprintf(stdout, "dhm optf.dhnifksp]); fprintf(stderr, "dvm optf.dvm[ksp]); fprintf(stdout, "dvii %6.3f; optf.dvm[ksp); fprintf(stderr, "dhn optf.dtm[ksp]); fprintf(stdeut, "dtun %63f~, optf.d[ksp]); fprintf(stderr, "dtc optf.dtc[ksp]); fprintf(stdout, "dhc %6.3f; optf.dhc[ksp); fprintf(stderr, "dhp optf.dhp[ksp]); fprintf(stdout, "dhp optf.dhp[ksp]); fprintf(stderr, "dvp optf.dvp[ksp]); fprintf(stdout, "dp %6.3f; optf.dvp[ksp]); fprintf(stderr, "dtp optf.dtp[ksp]); (ksp=O0; ksp <nsp; ksp++) S for ckk 0; ckk arsd[ksp].cns; ckk++) fprintf(stderr, "ksp ckk ",ksp, ckk); fprintf(stdout, "ksp ckk ",ksp, ckk); fprintf(stderr, "cnv ",arsd~ksp].cnv~ckkD); fprintf(stdout, "cnv arsd[ksplJ.cnv[ckk]); fprintf(stderr, !'dlt arsd[kspll.dlt[ckk]); fprintf(stdout, "dlt %6.3f; arsd[ksp].dlt[ckk]); fprintf(stderr, "wmn arsd[ksp].wmn[ckk]); fprintf(stdoutL "wmn arsd[ksp].wmn[ckk]); fprintf(stderr, "wxy arsd[ksp].wxy~ckk]l); fprintf(stdout, "wxy arsd[ksp].wxy[ckk]); fprintf(stderr, "wuv arsd[kspII.wuvfckk]); fprintf(stdout, "wuv arsd[ksp].wuv[ckk]); -103if initpr...d(&dd,&dn,&ki,&kj ,&nx,&ny,&ni,&nj ,&nih,&mv,nh,nv,kx,ky) exit( 1); if sdbrin..f(flmc, avrgm, imc: I, jmc kj, ni, nj) 0 exit(l); if sdbtimjf(flpc, avrgp, ipc ki, jpc kj, ni, nj) exit(1); if (itst =initprg.m(ki dn, kx, ky, nh, nv)

I

fprintf(stderr, "Subroutine initpr.... has failed, itst dn", itst); exit(itst); for 1; nv;j++) for (i 1; nh; hv~h0[k] i; vvOh0[kI =j; for (ksp ksp <nsp; ksp++) k= 1; for (j 1; nv; for i= 1; nh; arsd[ksp).uuvflk] ust; arsd[ksp].vvvf[0] vst; X. for (1lev lvs; 1ev >1Ist; 1ev--) sgdg =sgd; sgds SMTH sgd; sgmg =sgm; sgms SMTH sgmn; i~i fprintf(stdout,"The level sgrng sgdg 1ev, sgm*scl, sgd*scl); fprintf(stderr,"The level sgmg sgdg 1ev, sgm*scl, sgd*scl); scln 2.0 for ksp 0; ksp nsp; ksp++) arsdrksp].csq arsd[ksplj.c *arsdjjksp].c *scln; arsd[ksp].bsq arsd[ksp].b *arsd[ksp].b *scln sgmb;, arsd[ksp].psq arsd~kspl.p *arsd[ksp].p *scln; arsd[kspl.qsq arsdllksp].q *arsd[kspi~q *scln sgmnq; if rnxn continue; fprintf(stdout,'The derivatives of the images- fprintf(stderr,"The derivatives of the images: for (bn 0; bn bns; bn++) if (rdbtimjf(flmc, avrgm, optf.fgm, iroc ki, jmc kj, ni, nj, bri) !0) exit(l); fgg =optf fgm; ggg =optf.ggm[bn]; ggx optf.gxm[bri]; ggy =optf.gyrn[bn]; gxg optf.xgm[bn]; lnm -104gxx optf.xxm[bn]; gxy optf.xym[bn]; gsx optf.sxm[bn]; gsy optf~sym[bn]; setder(dd,dnjki~kj ,ni,nx,ny,fgg,ggg,ggx~ggy,gxg,gxx,gxy,gsx,gsy); fprintf(stdout,"fgm[%d], bn); fprintf(stderr,"fgmn[%dJ, bn).for (bn 0;t bn bns; bn++) if (rdbdtim(flpc, avrgp, optf.fgp, ipc ki, jpc kj, ni, nj, bn) =0) exit(l); fgg =optf.fgp; ggg =optf.ggp[bn]; ggx aptf.gxp[bn]; ggy =optf.gyplbn]; gxg optf.xgp[bn]; gxx =optf.xxp[bn]; gxy optf.xyp[bn]; gsx =optf.sxp~bn]; gsy optf.syp[bn]; setderdd,dn,i~j,ni,nx,ny,fgg,ggg,ggx,ggy,gxg,gxx,gxy,gsx,gsy); fprintf(stdout,"fgp[%d], bn); fprintf(stderr,"fgp[%d], fprintf(stdout, "have been set.\n"); Auk fprintf(stderr,"have been set.\n"); if corres(dddnkx,ky,nx,ny,mh,mv,nh, nvjrep,fep,mxn,mxu,mxl,rnxs) 1) exit(1); fprintf(stdout,'corres is completed.\n"); fprintf(stderr,"corres is completed.\n"); if prctof(&dd, &dn, &mnh, &mv, nh) 0) exit(l); ~:fprintf(stdout,'Trhe level %d is conipleted.\n", 1ev); Sfprintf(stderr,"The level %d is coinpleted.\n", 1ev); sgm sgin sqrt(msg); sgmb sgmb SGMB; e sgd sgd /sqrt(dsg); sgmq sgmq /SGMQ; rhx =rhx (sqrt(msg) sqrt(two)); :::for (lev =lstlev lev-) if prctof(&dd, &dn, &mnh, &mv, nh) 0) exit(l); fprintf(stdout,"The level %d is completed.\n", 1ev); fprintf(stderr,'Te level %d is completed.\n", 1ev); if wtdtfl(fldt, nh, nv, in'c, jrc, ip'c, jpc, kx, ky) =0) exit(0); initpr m(kp, dn, kx, ky, nh, nv) long kp, dn, kx, ky, nh, fly; float *gopt, *wrs *jj, *var, *work, *uujyf, *vvvf; long mars, ktnip, ksp, nsc, bn, i, j; ktmp kx dn; ktmp ktrnp dn; if bns 0 return(1); if (kx ktnp) fprintf(stderr,'Tbe kx should be divisible by dn); return(l); -105ktmp ky dn; ktmnp =ktmp cm; if (ky ktnp fprintf(stderr, tt rhe ky should be divisible by dn); retumn(l); ktnp (nh 1) /dn; ktmp =ktmp df+ 1; if (nh ktmp) fprintf(stderr,"The nh I should be divisible by dn); return( 1); ktrnp (nv 1) dn; ktmp =ktmp *dn +1; if (nv ktnp) fprintf(stderr,"The nv 1 should be divisible by dn); retun( 1); gopt (float malloc((unsigned)(16 nx bns sizeof(float))); if (!gopt) fprintf(stderr,"memory requested for array gopt failed, msp msp); *return(l); Ielse msp 16 nxy bns sizeof(float); fprintf(stderr,'gopt: msp msp); for (bn bn <bns; bnHoptf.ggm[bnl gopt (bns 0 bn) nxy; optf.gxm[bn] gopt (bns 1 bn) nxy; optf.gym[bn] gopt (bns *2 bn) nxy; optf.sxrnfbnl gopt (bns *3 bn) nxy; optf.xgm[bn] gopt (bns *4 bn) nxy; Am optf.xxm[bn] gopt (bns *5 bn) nxy; W optf.xym[bn] gopt (bns 6 6+ bn) nxy; optf.sym[bnl gopt (bns *7 bn) nxy; opfgg*n bs* n)*n optf.ggp[bn] gopt (bns 8 bn) nxy; optf.gxp[bn] gopt (brns 109 bn) *nxy; optf.gyp[bn] gopt (bns 11 1+ bn) *nxy; optf.sxp[bnl gopt (bns *12 bn) *nxy; optf.xgp[bn] gopt (bns *12 bn) *nxy; optf.xxp[bn] gopt (bns 13 bn) *nxy; optf.xyp[bn] gopt (bns *14 bn) *nxy; hv~hO (float rnalloc((unsigned)(nhv sizeof(float))); if !hvO&h0) fprintf(stderr," memory requested for array hv~h0 failed, msp dn", msp); return( 1); Ielse msp nhv sizeof(float); vvOhO (float malloc((unsigned)(nhv sizeof(float))); if !vv~hO) fprintf(stderr,"memory requested for array vvOhO filed, rnsp msp); remur( 1); }else msp nhv sizeof(float); uuvf (float I) malloc((unsigned)(nsp nhv sizeof(float))); if !uuvf) fprintf(stderr,"memory requested for array uuvf failed, msp msp); return(1); else msp nsp nhv sizeof(float); for ksp 0; ksp nsp; ksp++) arsd[ksp].uuvf uuvf +r ksp nhv; vvvf (float mafloc((unsigned)(nsp sizeof(float))); if !vvvf) fprintf(stderr,"memory requested for array vvvf failed, msp dn", msp); retun(l); Ielse msp asp sizeof(float); for ksp 0; ksp nsp; ksp++) arsd[ksplj.vvvf vvvf ksp; nsc =0; for (ksp ksp <nsp; ksp++) nsc nsc arsd[kspjl.cns; ktrnp (15 asp 4 nsc 17) *nhv *sizeof(float) (4 nsc 17) nhv *sizeof(long); ktmp nij sizeof(float) ktmp nij sizeof(float); uars (float malloc((unsigned)ktmp); if( uars) fprintf(stderr,"memory requested for array uars fadled, rasp dn", msp); retun( 1); else rasp ktmp; fprintf(stderr,"uars: rasp rasp); opff.fgm uars; optf.fgp uars; for (ksp ksp <nsp; ksp++) arsd[ksp] .u~vf uars; (nsp 0 ksp) nhv; arsd[ksp].udvf uars (nsp 1 ksp) nhv; arsd~kspl.utvf uars; (nsp *2 ksp) *nhv; arsd[ksp].duvf uars (nsp *3 ksp) *nhv; arsd[ksp].rucv 'uars (nsp *4 ksp) *nhy; arsd[ksp].pucv uars (nsp *5 I sp) *nhv; arsd~ksp].qucv uars (asp *6 ksp) *nhv; arsdlkspl.zucv uars (nsp 7 7+ ksp) *nhv; mrd[ksp].auul uars (nsp 8 ksp) rhv; arsd[ksp].auu2 uars (nsp *9 ksp) *nhv; arsd~ksp].g2sq uars; (asp *9 ksp) *nhv; arsd~ksp].auu3 uars; (asp *10 ksp) nhv; arsd[kspl.g3sq uars; (asp *10 ksp) nhv; -107arsi [ksp]. auu4 uars (nsp 1 i ksp) nhv; arsc[ksp].g4sq =uars +(nsp 11 ksp) nhv; uars (nsp 12 ksp) nhv; uars (nsp 12 ksp) nhv; arsd[ksp].ainv =uars (nsp *13 ksp) *nhv; arsd[ksp].dusq =uars (nsp *13 ksp) nhv; arsd[ksp].cuvl I uars (nsp *14 ksp) *nhv; wars =uars 15 *nsp *nhv; j=; for ksp 0; ksp 1 for i i arsd[ksp].cns; (rdkp.l~]=was+nc ~)nv arsdjjksp].w 1[i] wars (nsc 1 04 *nhv; arsd[ksp]. wl2rI i wars (nsc 2 1 nhv; arsd[kspl.w22[i] wars (nsc 3 2+j) *nhv; work wars 4* nsc nhv; wrkO0 work; woyk nhv; wrkOlI work; work nhv; wrk02 work; work nhv; wrk03 work; work nhv; wrkO4 work, work nhv; work; work nhv; wrk06 work; work nhv; ~:wrk07 work; work nhv; wrk08 work; work nhv; wrk09 work; work nhiv; :wrk 10 =work; work nhv; wrklIl work; work nhv; wrk 12 =work; work nhv; wrk 13 work; work nhv; wrkl14 work; work nhv; S wrkl15 work; work nhv; wrk 16 work; work nhv; mars =(lonig*) work; j =0; **for (ksp ksp <nsp; ksp+ for i 0; i arsd[ksp].cns; arsd[ksp].m IfIi] mars (nsc 0 j) *nhv; arsd[kspl.m12[i] mars (nsc I j) *nhv; arsd[kspl.m2 lfi] mars 4- (nsc 2 j) *nhv; arsd[ksp].m22[i] mars (nsc 3 j) *nhv; mars 4* ns nhv;, r~i'OhO mars; mars nhv; mvphp mars; mars nhv; mv~hp mars; mars nhv; mvmhp mars; mars nhv; -108invihO =mars; mars nhv; mvmhm =mars; mars nhv; mvOhm mars; mars nhv; mvphxn mars; mars nhv; mvphO mars; mars nhv; m 1m mars; mars nhv;, ml2m mars;imars nhv; m2 1 mn mars; mars nhv; m22m =mars; mars nhv; rnllp =mars; inars +--nhv; mIl2p =mas; mars +=nhv; in2lp =imars; mars nhv; m22P =rmn; vars =(float malloc((unsigned)(10 iip *sizeof(float))); if (!vars) fpririf~stderr,"inemory requested for array vars failed, misp rnsp); return(l); OA else insp 10 Pnsp sizeof(float); for ksp 0; ksp nsp; ksp++) (rdkp.~f=vas+np0+kp arsdllksp].vdvf vars nsp I 0+ ksp; arsdjlksp].vdvf= vars nsp 2 ksp; *ars-d[ksp.vvf vars +nsp32+ ksp; ars[kspl.vvfvars nsp 4 ksp; arsdrkspl.vcv =vars +nsp54+ ksp; arsd~kspj.pvcv vars nsp 6 ksp; arsd[kspl.zvcv vars nsp *7 ksp; arsd[ksp].bvvl vars nsp *8 ksp; arsd[ksplj.binv vars nsp 9 ksp; fprinrffstderr,"msp, nij nhv nxy msp, flu, nhv, nxy); return(0); sdbtimj'(fileln, averge, kid, kkj, ni, nj) char *filein; *::*float *averge; **.:long kki, kkj, ni, nj; unsigned char buffer[BUFEL]; float pxltln; long pixels, lines, bands, form, imgin; long apixel, aline, aband, arrtyp, anrtyp JDBYTE; imngin 0; if opnimg(&imgin, filein, 0, FALSE) I= SYSNRM return(l); getdef(&imgin, &pixels, &Mines, &bands, &form); nxltln pixels lines; if ki dk <0 dki =-kki; if (dki <kki +ni -pixels )dki kId+ a ,xels; if (dkj +kkj <0 likj= -kkj; if (dkj <kkj +nj -lines dkj kkj +nj -lines; if (form arrtyp) fprintf(stderr,"The input is the wrong data typefn"); -109return( 1); for aband aband bands; aband+-iaverge[aband] 0.0; for aline 0; aline lines; aline++) for aband 0; aband bands; aband+-irdline(&imgin, alixie, aband, buffer, BUFEL, arrtyp); for apixel 0; apixel pixels; apixel++) averge~aband] averge~aband3 buffer[apixel]; for (aband 0; aband bands; aband-Haverge[aband] averge[aband] pxltln; clsimg(&izngin); if bns 0) bns bands; else if bns !=bands return(1); return(0); *:~:rdbdtmj(filein, averge, fg, kki, kkj, ni, nj, bn) char *filein; float *averge, *fg; long kid, kkj, ni, nj, bn; unsigned char buffer[BUFEL]; long nirnd, pixels, fines, bands, line, band, imgin; long njmd, dpixel, dine, pixel, form, arrtyp, k; nimd =ni did; njmd nj dkj; 0arrtyp =IDBYTE; imgin =0; if opnimg(&imngin, filein,.0, FALSE) SYSNRM) return(1); getdef(&irngin, &pixels, &Mines, &bands, &form)-, if (form!=arrtyp) fprintf(stderr,"The input is the wrong data type~n"); return(l); if bands bns) fprintf(stderr,"The input is the wrong bands type\n"); return( 1); if bands BANDS) fprintf(stderr,"The input has too many bands~n"); retuin( 1); if pixe.Is BUFEL) -110fPrintf(std(-rr,"The buffer is not large enough\n"); return(1); if kki dki <0) fprintf(stdlemT"The did %d is not large enough, di); fPrintf(stderr,"the value of dkI should be k~d); return(l); if pixels ki nimd <0) fprintf(stderr,"The did %d is not large enough, the did); fprintf(stderr,"value of di should be ki ni pixels); return(l); if kkj dkj <0) fprintf(stderT,"The dkj %d is not large enough, dicj); fprintf(stdlerr,"the value of dkj should be -kkj); return(l); if (lines -kkj -njmd 0 fprintf(stderr,'The dkj is not large enough, the dkj); fprintf(stderr,"value of dkj should be kkj nj lines); return(l); for line kkj dkj; line klkj njmd; line++) for band 0; band bands; band++irdline(&imgin, line, band, buffer, BUFEL, arrtyp); if band =bn) 4W.: if line =kkj dkj k 0; for pixel ki di; pixel kki nimd; pixel++) fg[k] bufferf pixel] averge[band]; if dki>0) if ((pixel- k=ki+ did) 11 "pixel =kid nimd 1) for dpixel 0; dpixel did; dpixel++) fg[k fg[k if dkj 0) if (line =kkj dkj 11 line kkj njmd 1) for dline 0; dline dkj; dline++)

-III-

for pixel 0; pixel ni; pixel++fg[k] =fg[k nil; clsimg(&imgin); return(0); float *rtgvtr, *rtyvtr, *rggyfl., *rgxvLtr, *r.x trl, *rxvftr, *tggt, *tgyvt initpr..d(dd, dn, I, kj, nx, fly, ni, nj, nih, my, nh, fly, kx, ky) long *dd, *nx, *ni, *1221, nli, kx; long *dn, *kj, *ny, *nj, *pmy, nv, ky; double sqrtl, sqrt2, ceilO0, sqirQ; long ddv, kiv, nxv, niv, kgi 0, ksi 0, kgit2pl1, ksit2pl1, ktp; *fe long dnv, kjv, nyv, njv, kgj 0, ksj 0, kgj t2pl1, ksj t2pl1, kc; ddv sgmb thre scl :dnv sgmq sgmw for(k= 1;k<lvs~k++) 1 sgrn sgm sqrt(msg); sgnlg sgrn; sgms SMiTH sgm; sgmb =sgmb *SGMB; sgd =sgd sqrt(dsg); sgdg sgd; sgds SMiTH sgd; sgrnq sgmnq SGMQ; scl =sci sqrt(two); rhx rhx sqrt(two msg); sgmw =sgmnw *SGM4W; if (ddv dnv) dnv=2* ddv; sqrtl ktp =thre 2.0 ceil((double)(0.5 *sqrtl *sgmg)); if ktp kgi kgi ktp; ktp thre 2.0 ceil((double)(0.5 *sqrtl *sgdg)); if ktp kgj kgj ktp; ktp thre 3.0 ceil((double)(0.5 *sqrtl *sgms)); if ktp ksi ksi ktp; 0:00: ktp, thre 3.0 ceil((double)(0.5 *sqr~t1 sgds)); if ktp ksj ksj ktp; else ddv =dnv; sqrt2 sqrt(two); ktp =thre 2.0 ceil((double)(0.5 *sqrt2 *sgmg)); if ktp kgi kgi ktp; ktp thre 2.0 ceil((deuble)(0.5 sqt* sgdg)); if (ktp >kgj )kgj =ktp; ktp thre 3.0 ceil((double)(0.5 sqrt2 *sgrns)); if ktp ksi ksi ktp; ktp thre 3.0 cefl((double)(0.5 sqrt2 *sgds)); if ktp ksj ksj ktp; -112kgit2pl 2 *kgi+ 1; kgjt2pl 2 *kgj 1; ksit2pl1 2 *ksi+ 1; ksjt2pl= 2 *ksj +1; kIv ksi dnv; *mhjj (nh 1) dnv 1; kjv ksj dnv; *mv =(nv 1) /dnv 1; nxv =tau*(lb 2 *kx 1; *kI=kiv; *nx =nxv; nyv tau (nv 1) 2 ky 1; *kj kjv; *iiy nyv; niv rho *(nxv 1) 2 kciv 1; *dd ddv; *n niv; njv =rho (nyv 2 *kjv *dn =lnv; *nj =njv; flu fliv *njv 1; nhv =nh 1; nxy nxv *nyv 1; rtgvtr (float rnalloc((unsigned)(ksjt2p I sizeof(float))); if !rtgvtr) fprintf(stderr,"memrnoy requested for array rtgvtr failed, msp msp); return(l); I else rnsp ksjt2pl sizeof(float); rtyvur (float rnalloc((unsigned)(ksjt2pl sizeof (float))); if( rtyvtr) fprintf(stderr,"memory requested for arry rtyvtr failed, msp nisp); return(l); Ielse msp ksjt2pl sizeof(float); rggvtr (float rnalloc((unsigned)(ksit2p 1 sizeof(float))); if(!rggvtr) fprintf(stderr,"rnemory requested for array rggvtr failed, msp msp); Ielse msp ksit2pl sizeof(float); rgxvtr (float malloc((unsigned)(ksit2pl sizeof(float))); if (!rgxvtr) fprintf(stderr, "memory requested for array rgxvtr failed, msp dn", msp); return(l); Ielse msp ksit2p I sizeof(float); W rxgvtr (float malloc((unsigned)(kgit2plI sizeof(float))); if(!rxgvtr) fprintf(stderr,"memory requested for array rxgvtr failed, msp W dn", msp); return( 1); Ielse msp kgit2pl sizeof(float); rxxvtr (float malloc((unsigned)(kgit2pI sizeof(float))); if( rxxvrr fprintf(stderr,"memory requested for array rxxvtr failed, msp msp); return(l); Ielse rnsp kgit2pl sizeoffioat); tggvtr (float malloc((unsigned)(niv sizeof(float))); if( !ggvtr) fprintf(stderr,"rnemory requested for array tggvtr failed, msp msp); re turn(1); Ielse rnsp niv sizeof(float); tgyvtr (float malloc((unsigned)(niv sizeof(float))); -113if( tgyvtr) fpintf(stderr, "memnory requested for array tgyvtr failed, msp msp); retun( 1); else msp niv sizeof(float); fprintf(stderr, msp msp); return(O); setder(dd, dii, ki, ki ni, nx, ny, fgg, ggg, ggx, ggy, gxg, gxx, gxy, gsx, gsy) long dd, dn, Id, kj, ni, nx, fly; float *fgg, *ggg, *ggx, *ggy, *gxg, *gxx, *gxy, *gsx, *gsy; Subroutine **setder** computes the values of the following image functions: ggg, ggx ggy, gxg, gxx, gxy, ggs, gxr,, by correlating the values of the initial image function cgg with the values of the following two-dimensional arrays: rgg, rgx, rgy, rxg, rxx, rxy, rgs, rxs. The array, ggg, is related to the array, cgg, as follows: ggg[x nx y] corresponds to cgglix] ni j[y]j where: i[xI=ki+rho*x, and j[y] =kj +rho for x=O, ,nx-l1, and 1. /~The arrays ggx, ggy, gxg, gxx, gxy, ggs, gxs, are related to thearray cgg in the same way as the array ggg is related to the array cgg. The variables are defined as follows.

Images (each represented as a 1-dimensional array, in a row major order): cggi +ni jI, ni j nj is the initial image value at the pixel (kid di i) on the line (kkj -dkj ggg[x +nx x nx y is I correlated image value at pixel i[x] (I rho x) on the I* line j [y]j (kj rho y) in the i, j coordinate system; 1* the correlation function rgg is defined to be a modified I* Gaussian kl-nction. 1* gxg[x +nx x=0, nx- 1, 1, is corkelated image value at pixel (i rho x) on the line j (kj rho y) in the i, j coordinate system; the correlation function rxg is defined to be a modified first-order partial with respect to variable x derivative of the Gaussian function.

ggx[x +nx x nx 1, y ny is the value of the first-order partial with respect to the variable x derivative of the function, ggg, where the system 1* of coordinates relates to the system of coordinates -114- 1* (ij) as: ifx] (ki rho j~y] (kj rho 1*ggy[x yI, x=O0, nx 1, y ny 1, 1* is the value of the first-order partial with respect to the 1* variable y derivative of the image, ggg, where the system of /~coordinates relates to the system of coordinates (ij) i[x] rho j[y (kj +rho 1* gxx[x +nx x=0, nx y=O0, ,ny -l 1, f~is the value of the first-order partial with respect to the variable::1 derivative of the function, gxg, where the system /~of coordinates relates to the system of coordinates 1* (ij) as: i[x] (ki+ rho j~y] (kj rho 1* gxy[x +nx x=O0, nx- 1, y ny /~is the value of the first-order partial with respect to the variable y derivative of the image, gxg, where the system of /~coordinates relates to the system of coordinates (ij) as: i[x] =(ki +rho j[y] =(kj +rho 1* ggs[x +nx* x nx-l1, y=O0, ny 1, 1* is the value of the smoothed vertion of the function, ggg, at the pixel i[x] (ki+rho on the line j[y] =(kj +rho gxs[x +nx* x nx-1, y=0, ,ny -l 1, is the value of the smoothed vertion of the function, gxg, at the pixel i[x] (I+rho on the line j[y] =(kj +rho Resolution pyramid: (dn[level], level ,..,levels), where: dn[l] 1, and for level levels dn~level] is defined as follows: 1* dn[level] =2 ddnlevel-l], if dn~level-1] dd[level-1] 1* dn [level] dd~level-lI], if dn(level-l]J= 2 *dd~level-l]; /*1 1* (dd~level], level levels), where: dd[1] 1, and 1* for level levels ddplevel] is defined as follows: dd[levelI dn[level-l], if dn[level-1] 2 dd[level-l] 1* ddlevel]j 2 dnolevel-1], if dnplevel-1] dd[level-l]; /~For each level 1, ,levels, the Gaussian parameters: kg1*vl =*1eci~gm)d~ee] 1* kgi[level] thre *ceil(sgmg) *dn[level], 1* ksj[level] thre ceil(sgdg) *dn[level], 1* ksi[level] thre ceil(sgms) *dn[level], ki ksi~levels]; kj ksj[levels); long k, nitdd; register float *fg, *tg, *ty, *rg, *Iry; float *gO, *gl, *g2, *g3, *g4, *g5, *sO, *sl; -115- L oat *rtg, *fly, *rgg, *rgx, *rxg, *rx *rsg, long i, ii, kgi, ksi, ngi, znrk, ifgidd, ifsidd, kgit2pl, ksit2pl; long j,jj, kgj, ksj, nsi, igg, ifgjdd, ifsjdd, kgjt2p 1, ksj t2p 1; double omneg, rtmp, dtmp, tnipxs, sqrtl, sqrto, sgmgi, sgmsi; double omsq, tmsq, dscl, trnpys, sqrt2, ceilO), sgdgj, sgdsj; nitdd =ni *dd, if dd =dn) sqrtl kgi dire 2.0 ceil((double)(0.5 sqrtl sgmg)); kgj thre 2.0 ceil((double)(0.5 sqrtl sgdg)); ksi dire 3.0 ceil((double)(0.5 sqrtl sgms)); ksj thre, 3.0 ceil((double)(0.5 sqrtl sgds)); else sqrt2 sqrt(two); kgi thre *2.0 ceil((double)(0.5 sqxtA2 sgmg)); kgj thre *2.0 ceil((double)(0.5 sqrt2 sgdg)); ksi thre *3.0 ceil((double)(0.5 sqrt2 sgms)); ksj thre *3.0 ceil((doubl6)(0.5 sqrt2 sgds)); orneg 3.0; ornsq omeg onieg; dtmnp rho tau; dscl dd /dtmp; :sgmgi sgmg *scl /dtnip; sgdgj sgdg *scl /dtmp; sgmsi sgms *sc dtnp; *:sgds sgds *scl dhp; *::rtg =rtgvtr-,k =0; rty rtyvtr, rtmnp 0.0; for (j=-kgj; kgj; :trnpys =j *dscl /sgdgj; tmsq tmpys tmpys; dtmp -0.5 tnisq (1.0 omsq *rs) rtg =exp(dtmp); rtmp rtmp rtg[k]; trnsq =1.0 omsq ttnsq; rty[k] =rtg[k] tmpys (sgdgj *tinsq *tmsq); rtnp sqrt((double)isc) rtmp; k =0; for (j =kgj; j kgj; rtI ]=rgk tp rtg[k] rtgllk] rtmp; rtmp 0.0; k 0; rgg rggvtr, rxg rxgvtr, rgx rgxvtr, rxx rxxvtr, for i =-kgi; i kgi; tnipxs i dscl sgnmgi; tmsq tmpxs tmpxcs; dtnp -0.5 tmq (1.0 omsq *trnsq); -116rgg[kI exp(dtrnp); rtrnp rtrnp rggllk]; tmnsq (1.0 omsq tmsq) (1.0 omsq tmsq); rgx[k] rgg[kI tmpxs (sgmgi *tmsq); rxg[k] rgglk] tmpxs sgmgi; rxx[k] rgg[k] (trnpxs tmnpxs rmsq (sgmgi *sgrngi); k++ rtmp sqrt((dotible)isc) rtmp; k =0; for i =-kgi; i kgi; i-Hrgg[k] rgg[kI r t mnp; rgx[k] rgx[k] *rtmnp; rxg~k] rxg[k] *rtnip: rxx[k] rxx[k] *rhrlp; if dn dd xnrk 0; else mrk 1; kgit2pI 2 *kgi kgjt2p I= 2 *kgj +1; fgjdd =ni *(kj -kgj *dd) +kI-kgi dd; igg=0; ngi =rho (nx 1) 2* kgi 1; for (j 0; j ny; j dd) ifgidd =ifgjdd; t9g=tggvtr y tgyvtr, for (i i<ngi; dd) fg =fgg +ifgidd; **tgo= rgl= tg; *ty =00; ry ty for (j 0;jj kgjt2pl1;jj++) *tg *ty fg +=nitdd; ifgidd dd; tg-H-; ty++; ifgjdd ni *rho *dd; tg =tggvtr, ty tgyvtr, gO ggg igg; gl ggx igg; g2 ggy igg; g3 =gxg +igg; g4 =gxx +igg; g5 =gxy +igg; for i 0; i nx; i dd) if mrk 0) rgg rggvtr, rgx rgxvtr, rxg rxgvtr. lxx rxxvtr, *gO 0.0; *gl 0.0; *g2 0.0; *g3 0.0; *g4 0.0; *g5 0.0; for ii 0; ii kgit2p 1; ii++) *go *tg *rgg; *gl *tg *rgx; *g2 *ty *rgg; *g *tg rg *g4 *tg rx -117- *ty *rxg; rgg++; rgx++; tg++; rxg++; rxx++; ty++; tg kgit2pl1; ty kgit2p 1; tg rho; ty rho; rnrk =-mrk; gO dd; g I dd; g2 dd; g3 dd; g4 dd; g5 dd; igg nx dd; rtrnp 0.0; k 0; for (j =ksj; ksj; tmpys j dscl sgdsj; dtmp tmnpys tmnpys; Srtg[k] exp(dtmp); rtmp rump +rtg[k];, W rty[k] =rtgfk] txz;'ys sgdsj; k)+ rtmp sqrt((double)isc) rtmp; k =0; j=ksj; ksj; rtg[k] =rtg[k] *rtmp; Itk]= rty[k] rtmnp; rsg rggvtr, rsx rgxvtr, rrnpxs i*dscl /sgnsi; dtmnp trnpxs tmpxs; r sg[k] =exp(dtmp); rtxnp rtmp rsg[k]; sx[k] rsg[kI tmpxs sgmsi; sqrt((double)isc), tr k 0 *fcr( i=ksi; ksi; i++ rsg[k] rsgfk] r t tnp; rsx[k] rsx[k] *mrtnp; k-H-; if &I dd mrk 0; else nirk 1; ksit2pl =2 ksi+ 1; ksjt2pl =2 *ksj +1; ifsjdd ni (kj ksj dd) ki ksi dd; igg =0; nsi =rho *(nx 1) +2 *ksi 1; for j 0; j ny; j dd) ifsidd ifsjdd; t9 tggvtr, ty tgyvtr, for( i 0; i <nsi; i dd) m -118fg =fgg+ ifsidd; *tg 0.0; rg rtg; *ty=O.0;ry ty; for~j jj ksjt2pl1;jj++) *tg rg *ty *fg *ry+ nitdd; ifsidd dcl; tg+-i; ty++; ifsjdd ni rho *dd; sO gsx igg; tg =tggvtr, s I gsy +igg; ty tgyvt, for (i 0; i nx; i dd) if mrk 0 rsg rggvtr- *s0 0.0; rsx =rgxvtr, *slI 0.0; for ii 0; ii ksit2p 1; ii++) *s rx 9. *s1 sg+ tg -=ksit2p 1; ty ksit2pl1; mrk rrk; sO dd; tg +=rho; slI+= dd; ty+= rho; igg fix dd; .ewres(dd, dn, kx, ky, fix, ny, nih, my, nh, fly, reps, feps, rnxn, mxii, mxl, rnxs) '!rMng dd, dn, kx, ky, nx, ny, m~h, my, nh, nv, mxn, mxu, rnxl, mixs; *:tlugt reps, feps; register long k, kOhO; long ksc, nsc, itn, itu, its, itst; float rkf, ff0, stp STEP, stl, std, fepsg, tmpuu, tmpvv, temp; std STPD; if itst setinl1(dd, dn, mh, my, nh)) =0 tjrintf(stderr,"Subroutine corres: setin 1") fprintf(stderr,"returnes value equal to itst); return(itst); if (rnxn<lI return(0); for ksp 0; ksp nsp; ksp++) for (k =kv~h0; k; k--

I

kOhO mv~h0[k]; -119arsd[ks 1 4.u~vf[kOhO] =arsd[ksp].uuvf[kOhO]; arsd[ksp].v~vf[OI arsd[ksp].vvvf[O]; for itn 0; itn mxn; itn++) for (.ksp 0; ksp nsp; ksp+ if ((itst setdr2(dd. dn, kx, ky, nx, ny) 1= 0) fpr~ntf(stderr,"Subroutine corres: setdr2 fprintf(stderr,"returnes value equal to UdM", itst); return(itst); rf0 0.0; for ksp ksp <nsp; ksp++) if itst setmtx(dd, dn, nh, nv, &rfO) 0) fprintf(stderr,"Subroutine corres: setmtx fprintf(stderr,"returnes value equal to itst); return(itst); foi( k =kv~h0; k; karsd[ksp].udvf[xnv~h0[k]] 0.0; *::arsd[ksp].vdvf[0] 0.0; :.fepsg =feps *rf0; :fpintf(stderr, "ri' rk ,rf0 (float) (nsp IvOhO)); fiprintf(stdout, "rf rk ,rf0 (float) (nsp *kvOh0)); for its 0; its mxs; its++iconjgr(&rkf, fepsg, rnxl); fprintf(stderr, rkf/ (float) (nsp *kvOh0)); fprintf(stdout, rkf/ (float) (nsp *kv~h0)); *:kmrpuu 0.V: tmpvv 0.0; for ksp ksp nsp; ksp++) for kkvOh0; k; kOhO =mv~h0[k]; arsd[ksp.uuvf[kOhO] arsd[ksp.uuvf[kOhOI stp arsd[ksp].udvf[kOhO]; lnpuu arsdj[ksp].uuvflkOh0]; arsd[kspl.vvvf[0] arsd[ksp].vvvf[0] sip arsd[ksp].vdvf[0]; tnipvv arsd[kv;p].vvvf[0I; temp nsp kv~h0; tmpuu tmnpuu temp; temp nsp; tmpvv trnpvv temp; -120fprinwf(stderr, tmpuu, tmpvv); fprintf(stdout," tmpuu, trnpvv); return(0); prctof(dd, din, nih, mv, nh) long *dd *dxl, *mh, *mv, nh; register long k, k1 1, k12, k21, k22; long nhpl1, ki1, i, j, dntnh, ddtnh, ddtnhp 1; if (*dn (*dd 1) return(O); nhpI nh+ 1; dntnh =*dn nh; if *cln *dd) *dd *d ddtnhpl I *dd nhpl; k =kl; kl dntnh; for( i 1 *flh; k11I=k -ddtnhp1; kl2=k11+ *dn; k22 k +ddtnhpl; k21= k22 -*dn; *:~for ksp 0; ksp nsp; ksp++) arsd[ksp].!'nvvf[k] (arsd[ksp].uuvflkl 1] arsd[ksp].uuvf[kl 2] arsd[ksp].uuvf[k21] arsd[ksp].uuvflk22]) *0.25; k. *dn; *Itturn .:l*dn=*dd +*dd) *Udtnh *dd nh; ki *dd+ 1; k =kl; kl dntnh; kIlI *dd; k22= k +*dd; for ksp 0; ksp nsp; ksp++) arsdfksp].uuvflk] (arsd[ksp].uu%. 1k 1 1] arsd~ksp].uuvflk22]) k *dfl; for (j j j+4 k =kl; ki dntnh; -121for( i= 1; j <*mhJ; k12 =k ddtnh; ki 1I k *dd; k2l k +ddtnh; k22 k +*dd; for ksp 0; ksp nsp; ksp-H-) arsdlksp].uuvfllk] (arsd[ksp].uuvflkl 1] arsd[ksp].uuvflkl2l arsd[ksp].uuvflk2l] arsd~ksp].uuvf[k221) *0,25; k *dn; k kl; for( i 1 i *mhj; ij+) kl *dd; k2=k +*dd; arsd[ksp].uuvfljk] (arsd[ksp].uuvflkll] arsd[ksp].uuvffk22]) k *dn; kl =*dd *nh +1; (j 1;j k =k1; klI dnmnh; k 12 k ddtnh; k21 k ddmnh; for (ksp ksp <nsp; ksp++) arsd[ksp] .uuvf[k] (arsd[ksp].uuvflkl2] arsciksp].uuvf[k2l]) for(i=2; k12 k ddtnh; ki 1I k *d k21 k ddtph; k22 =k *dd; for ksp 0; ksp nsp; ksp++) arsd[ksp].uuvffk] (arsd[ksp].uuvf[kl 1) arsd[kspj.uuvflk1l2] ars-dtksp].uuvffk2l] arsd[ksp].uuvftk22]) 0.25; k *dfl; k 12= k -ddtnh; k2l k +ddtnh; for (ksp ksp <nsp; ksp++) arsdlkspl.uuvflk] (arsdfksp].uuvflkl2l arsd[ksp].uuvffk2l)) -122- *dn *dd; 2 (*mv J) 1; return(0); return~l); wtdtfl~fldt, nh, nv, kim, kjm,, ip, kjp, kx, ky) char *fldt; long nh, nv, im, kjm, ip, kjp, kx, ky; float tenmp, buffer[BUFEL]; int wrtn; static int fd; long k, i, ix, jy, bn, ksp, nbuffer, fd open(fldt, QRDWR I 0_GREAT, 02644); Of(fd<0) fprintf(stderr, "Can't open file %Anz", fidt); retum( 1); :::Jemp nh; :irtn writed &mp sizeof(float) 4 pint sterr "write has failed at pon Wt close(fd); *return(l); *ftopnv; :::-wrtn =write(fd, &temp, sizeof(float)), wrtn 1= sizeof(float)) S. fprintf(stderr, "write has failed at point clo~e(fd); return(l); tnp =kin rho x 0 w-tn =write(fd, &temp, sizeof(float)); if wrtn. I= sizeof(float))

I

fprintf(stderr, "write has failed at point 4\n) close(fd); return(l); temp =kjm rho* wrtn =write(fd, &temp, sizeof(float)); if wrtn si-7eof(flzat))

I

fprintf(stderr, "write has failed at point close(fd); return(1); temp =kip rho kx; wi-t= S &temp, sizeof(float)); if w izeof(floa*,)) fprintf(stderr,vwrite ias failed at point close(fd); ret urn(1); teml kjp rho ky; wrm write(fd,&temp,sizeof(float)); if wi-tn !=sizeof(float))

I

fpi-intf(stderr," write has failed at point close(fd); return(l); temp =rho; wi-tn =wi-ite(fd,&temp,sizeof(float)); A*f wi-tn sizeof(float)) fpiintf(stderr," write has failed at point close(fd); retun( 1); :,I.mp =mtu; '.whn =wi-ite(fd,&temp,sizeof(float)); wirn !=sizeof(float)) f~prirnff(stderr,"write has failed at point IOMn"); retun( 1); n nsp; wmt vi-ite(fd,&temp,sizeof(fI-oat)); if wi-tn sizeof(float)) %'!'p-intf(stderr,' write has failed at point close(fd)X ieu l) ksp=; ksp <nsp; ksp++) temp =optf.dtinlksp] optf.dtc[ks-J; wrtn =write(fd,&temp,sizeof(float));, if wrtn sizeof(float)) fprintf(stderr,"wi-ite has failed at the point 18: fprintf(stderr,"temp wi-tn ",temp, wi-tn); fprintf(stderr,"ksp ksp); close(fd); return(1); temp =optf.dtp[ksp] optf.dtc[ksp]; wi-tn =write(fd,&temp,sizeof(float)); if wi-tn sizeof(float)) -124fprinrf(stderr, "write has failed at point 19.\n" t close(fd); return(1); k=1; nbuffer =4 *nh; if (nbuffer BUFEL fPr'ntf(stderr,'The buffer is not large enough~n' t close(fd); return(1); fo~r(jy 0; jy nv; jy++) i =0; for (ix 0; ix nl; ix++) bufri v(Ok;i+ buffer[i] whvO[k]; buffer[i] vvrh[k];uv~];i+ buffer[i] arsd[ksp].uuvf[k]; bkfri]=asdkp.vv;];i+ :~wrtn write(fd, buffer, nbuffer sizeof(float)); if~ wrtn nbuffer sizeof(float)) fprintf(stderr,"write has failed at point 23 An"); return(l); jbse(fd); *--ielurn(0); *ai2(dd, dn, kx, ky, nx, ny) long dd, dn, kx, ky, nx, fly; Oppble temp, tnipl, tmp3, expo; register long -w *fidat *hcp, *vcp, *ggpm, *gxm, *gym, *xgm, *xym, *sxm, *sym, rhgcn, scd; float *hcm, *vcm, *ggp, *gxp, *gyp, *xgp, *xxp, *xyp, *sxp, *syp. rhxcn, scn; float *dxlp, *dylp, *dy2p, *duul, *auul, *cuvl, gug gusq, gvsq, shx; float *Cixlm, *dylm, *dx~m, *dy2m, *dwl, *bwl, *jycn, gvgt, gugv, gtsq, shy; float temp 1, tmpgm, tempu, tmpgt, tmnpgx, tmpxx, tmpsx, *tmpg2, *tmpg4, tempw; float temp2, tlmpgp, tempv, tmpxt, tmpgy, tmpxy, tmpsy, *tmpg3, *tmpg5, *dusq; Iloat *uuvf, dxlpk, dx2pk, dylpk, dy2pk, dtp, dtc, dtrn; float *vvvf, dxlmk, dx2rnk, dylmk, dy2mk; lung kI Ip, kl2p, k2lp, k22p, bn; long ki lm, kl2m, k2lm, k22m, it; hcp wrk~O; vcp wrk~l; hcm wrk02; vcm wrk03; dxlIp =wrk04; dy Ip wrk5; dx2p wrk06; dy2p vvrk07; dxlr wrkO8; dyim wrkO9; dx2m AIrk; dy2m wrkl 1; lu-npg2 =wrkl2; tmnpg3 wrkl3; tznpg4 wvrkl4; tmpg5 dusq arsd[kspl.dusq; bcn wrkl6; uuvf arsd[ksp] .uuvf; vvvf arsd[ksp].vvvf, -125scd 2.0 dd dd; sed 1.0 /scd; scn dn *dn; scn 1.0 /scn; for( k =kvOh0; k;

I

dusq[mvOh0[kll 0.0; for k.=kvphp; k; tempu =uuvffmvmhmnfk]] uuvffmvphp[k]l; tempw =tempu tempui scd; dusq[mvmihm[kjj tempw; dusq[mvphp[k]] tempw; for( k =kv~hp; k; ktempu =uuvf[mv~hm[k]] uuvf[mv~hp[k]]; tempw =tenmpu tempu scn; *dusq[mv~hm~k]] tempw; Wdusq[znvOhp[k]j teinpw; for( k =kvmihp;k; k-- .tempu =uuvflrnvphm[k]] uuvf[mvmhp[k]]; empw =tempu ternpu scd; dusq[mvphm[k]] texnpw; .:dusq[mvmhpfk]] tempw; fr k =kvnh0; k; tempu =uvlvf~mvphO[k]] uuvf[mvmhO[k]]; -tempw tempu *tempu *scn; u&A[mvph0[k]] tempw; .*:dus[mvnmh0[k]J tempw; a;u I arsd[ksp].auu 1; Itv I arsd[ksp].cuvl; bvv I arsd[ksp].bvvl; :4muI arsd[ksp].duvf; jyl= arsd[ksp].dvvf; *fdr( k =kv~h0; k; kin mvOhO[kl; auulfml 0.0; cuvlI[m] 0.0; duul[m] 0.0; tinpg2f m] 0.0; tinpg3f m] 0.0; tmnpg4[mjj 0.0; 0.0; 0.0; 0.0; dtp optf.dtpfksp]; dtc optf.dtc[ksp]; dtin optf.dtm[kspl; -126for( k =kv~h0; k; kin mv~h0[k]; bcn[m] -hcp[m] =hv~h0[m] optf.dhp[ksp] (drp dtc) uuvf~m]; hcin[m] hv~h0[m] cptf.dhm[ksp] (dtm dtc) uuvf[rn]; vcp[mJ vOh0[m] optf.dvp[ksp] (dtp dtc) vvvf[0); vcm[m] vv~h0[i] optf.dvm[ksp] (dtm dtc) wvrt[0]; it =setin2(bcn,hcp,vcpmxll1p,ml2p,m2 lpm22p,dx ip, dy lp,dx2p,dy2p,dd,dn,kx,ky,nx,ny); if (it 0 0) return(it); it =setin2(bcn,hcinvcinil lrni,ml 2m,n2 lm,in22m,dx 1 i, dy 1 i,dx2in,dy2m,dd,dn,kx,ky,nx,ny); if (it! return(it); for bn 0; bn bns; bn++) ggm optf.ggin[bn]; ggp optf.ggp[bn]; Ab.gxin optf.gxm[bn]; gxp optf.gxplbn]; gym optf.gyin[bn]; gyp optf.gyp[bn]; xgm optf.xgxn[bn]; xgp optf.xgp[bn]; xxin optf.xxm[bn]; xxp optf.xxp[bn]; xyin optf.xyin[bn]; xyp optf.xypllbn]; sxin optf.sxxn[bn]; sxp optf.sxp[bn]; syin optf.sym[bnl; syp optf.syp[bn]; *for( k =kv~h0; k; km mv~h0[k]; **kl lp =ml lp[m]; dxlpk =dxlp[m]; I 2p In2p[in]; dx2pk =dx2p[m]; k2lp in2lp[in]; dylpk dylp[m]; k22p n2.2p[in]; dy2pk dy2p[m]; kl lmn=iml lwm]; dxlmk =dxlm[m]; Il2m inl2in[m]; dx2ink dx2m[m]; k2lin in2lin[m]; dylmk dylm[m]; k22m mn22-1n] dy2mk dy2m[m]; .:temp 1 dx2pk ggp[kll1p] dxlpk ggp[kl2p]; W temp2 dx2pk ggp[k2lp] dxlpk ggp[k22p1; txnpgp =dy2pk teinpi dylpk temp2; tempi dx2mri ggin[kllm] dxlmk ggin[kl2m]; *:~:teinp2 dx2ink ggm[k2lm] dxlmk ggin[k22in]; tmpgm =dy2nk temp dylnk teinp2; tnipgt tmnpgp tinpgin; teinpl dx2pk gxpfkl lp] dx lpk gxp[kl2p]; temp2 =dx2pk gxp(Ic2lp] clxlpk gxpllk22p]; tmpgp =dy2pk teinpi dylpk temp2; tempi dx2rnk gxm[kl urn] dxlmk gxm[kl2nj; teinp2 =dx2ink gxm[k2lin] dxlmnk gxrn[k22m]; tmpgrn dy2ink temp 1 dylmnk temp2; tmpgx dtp tinpgp dtrn tnipgm; temp I dx2pk gyp[kll1p] dxlpk gyp[kl2p]; emp2 dx2pk gyp[kalp] dxlpk gyp[k22p]; trnpgp dy2pk temp 1 dylpk teinp2; teinpi dx2rnk gyni~kllIm] dxlink gym[kl2in]; teinp2 dx2ink gym[k2lin] dxlmk gym[k22m]; tmpgm: dy2ink temp 1 dylink temp2; -127tmpgy dtp tmpgp dtm tmpgm; templ dx2pk *xgp[kl ip] dxlpk xgp[kl2pI; temp2 dx2pk *xgp[k2lp] dxlpk xgp[k22p]; tmpgp =dy2pk templ1+ dylpk temp2; tempi =dx2mk xgm[kllm] +dxlmk *xgm[kl2m]; temp2 dx2mk xgm[k2lm] dxlmk xgm[k22m]; tmpgmn dy2mk temp I dy Imk temp2; trppxt tmnpgp tmpgm; tempi dx2plc xxp[kl 1p] dxlpk xxpflkl2p]; temp2 =dx2pk xxp[k21p] dxlpk xxp[k22p]; tmpgp =dy2pk temp 1 dylpk temp2; temp 1 dx2mk *xxm[klIim] dxlrnk xxm[kl2m]; temp2 =dx2mk *xxmnlk2lm] dxlmk xxm[k22m]; trnpgm dy2mk temp I dylImk temp2; tmpxx dtp tmnpgp dtm tmpgm; tempi dx2pk xyp[kl 1p] dxlpk xyp[kl2p]; temp2 dx2pk xypfjk2lp] dxlpk xyp[k22p]; tmpgp dy2pk tempi dylpk temp2; temp 1 dx2mk *xym[kI lm]l dxlmk xym[k12m]; temp2 dx2mk *xym[k2lm] dxlmk xym[k22m]; tmpgm dy2mk temp 1 dylmk temp2; tmpxy =dtp tmpgp +dtm trnpgm; rhgcn rhg bcn[m]; rhxcn rhx bcn[m]; ****gtsq =rhgcn tmpgt *tmpg rhxcn tmpxt rnpxt; gugt rhgcn *tmpgx *tmpgt rhxcn tpxx tmpxt; gvgt rhgcn *tmpgy *tmnpgt rhxcri =pxy tmpxt; gus. rhgcn tmnpgx tmpgx +rhxcn tnpxx tmpxx; gugv =rhgcn *tmpgx *tmpgy rhxcn tmnpxx *tmnpxy; gvsq rhgcn *tmpgy *tmpgy rhxcn trnpxy *tmpxy; temp =(arsd~ksp].psq arsd[ksp].qsq arsdfkspil.dusq[m]) *gtsq; tmpl =arsd[ksp].r exp((double)(-I.0 temp)); *atmp3 =arsd[ksp).r exp((double)(-3.0 temp)); :duu I[m] -=gugt *tmp1; dvv1 -=gvgt *tmp 1; auulI[m] gusq tmp3; cuvl1[m] gugv tmp3; gs m3 tempi dx2pk sxp[kl lp] dxlpk sxp[kl2p]; temp2 dx2pk sxp[k2lp] dxlpk sxp[k22p]; ~:tmpgp =dy2pk temp dylpk temp2; temp 1 dx2mk sxm[kl Intl dxlmk sxmllkl2m]; temp2 dx2mk sx-m[k2lm] dxlrnk sxmfk22m]; tmpgm dy2mk tempi dylmk temp2; tmpsx dtp tmpgp dtm tmpgm; temp 1 dx2pk syp[kl ip] dxlpk syp[kl2p]; temp2 dx2pk syp[k2lp] dxlpk syp[k22p i; !rnpgp dy2pk tempi dylpk temp2; tempi dx2mk symlklm] dxlmk sym[kl2m]; temp2 =dx2mk sym~k2lm] dxlmk sym[k22m]; tmipgm dy2mk temp 1 dy Imk temp2; trnpsy =dtp tmpgp +dtmn* txnpgm; temp tmpsx tunpsy; tmpg2[m] 0.5 temp temp; temp tmpsx; timpg3[m] temp temp; temp tinpsx tmpsy; tmpg4[m] 0.5 temp temp; temp tmpsy- tmpg5[m] temp temp; -128for( k =kvphp; k; arsd[ksp].g2sq[mvph tmpg2[mvphp[k]] tmnpg2[mvmhm[k]I; for (k =kvOhp; k; arsdfksp] .g3sq[mvOhp[k]] tmnpg3[mvOhp[k]] tmpg3[mvOhmfk]I; for( k =kvmhp; k; karsd[ksp].g4sq[mvmhp[k]] tmpg4[mvmhp[k]] tmnpg4[mvphm[k]]; for k kvmhO; k; k-) arsd[ksp] .g5sq[mvmhO[k]] tmnpg5[mvmhO[k]] .**Vtin2(bcn, hc, vc, mal1, m12, m21, m22, dxl, dyl, dx2, dy2, dd, dn, kx, ky, nx, ny) *..Ioh'g *ml 1, *m12, *rn2l, *n22, dd, dn, kx, ky, nx, fly; .:.Toat *bcn, *hc, *vc *dxl, *dyl, *dx2, *dy2; ThIis subroutine perform the following functions. 1. The vector field I k =,.,kvOhO) M- is transformed into the vector field (xc I1k 1, kvOhO)* using the relations xc[k] =kx +tau *(hc k =1,..kvOhO, 1*yc[k] =ky +tau (vc[k] k =1,..kvOhO. 1* 2. Integer valued fields 1* (ml m12[k],nm.21[k], m22[kJ I k= 1..kvOhO)* I' are determined such that the metrix 1* m21 m.22[kI 1* mll~k] m12[k] formns a grid cell for every k 1, kvOhO. 1* 3. Real valued -fields -129- 1*x I dylI[k],dx2[k, dy2[k] I k 1, kvOhQ)* are determined such that the function value is a bilinear interpolation of its values at the grid points 1* f(mllI[k]), f(m1 f(rn21.[k]), f(m22[k]) 1* defined by the relations 1* f(xc[k],yc[k]) dy2[k] (Wx[k] f(ml 1 dx I[k] f(m 1* dy I[k] (dx2[k] f(mp2I dx double ficoro; float thre, hri, xm-int, xmaxt, ynunt, ymaxt; long mx, my, mxx, myy, nxtdm, nxtdp, nxtdn; float xmin, xmax, ymin, ymax, *xc, *yc, xcr, ycr; *xc =dxl; yc =dyl; thre =tau *dn; thri 1.0 /thre; xmin xmint xmin thre; xmax =nx kx 1; xmaxt xmax thre; ymin =-Icy; ymint ymnin thre; *.ymax ny -ky yaxt =ymax -thre; ::(register long k, kOhO; .:..'register float *rxc, *rhc, *ryc, *rvc; rc= xc; rhc hc; ryc yc; rvc vc; ror kkvoho; k; kkOhO =mvOh0[kIJ; i.:rxc[kOhO] tau *(rhc[kOhO] *::ryc[kOhO] tau *(rvc[kOhO] if (rxctlk~h0]<= xmin) bcn[kOhO] 0.0; :else if rxc[kOhO] xmint) bcn[kOhO] bcn[kOhOI (rxc[kOhO] xmnin) *hr -bcn[kOhOI (rxc[kOhO] xnuin) hi (rxc[k~h0]>= xmax) *::bcn[kOhO] 0.0; else if rxc[kOhOI xmaxt) bcn[k~hO] bcn~k~h0] (xxnax rxc[kOhO]) tMr bcn[kOhO] (xmax rxc[kOhO]) hi if ryc [kOhO] ymin) bcn[kOhO] 0.0; else if ryc[kOhOl ymint) bcn[kOhO] bcn[kOhO] (ryc[kOhO] ymin) tlu-j bcn[kOhO] (ryc[kOhO] ymin) thri; if ryc [kOhO] ymax) bcn[kOhO] 0.0; else if ryc[kOhO] ymaxt) bcn[kOhO] bcn[kOhO] (ymax ryc[kOhO]) thri -bcn[kOhO] (ymax ryc[kOhO]) thri; if dii dd)

I

-130nxtdn =nx clw xmrin Ix; xrnax nx kx 1; ymrin=-ky; yrax =ny -ky 1; (register long k, kOhO, *rml 1, *rm12, *rmr2l, *rm.22; register float *rc *rdx 1, *rdx2; register float *ryc, *rdy 1, *rdy2; rml I ml rml2 m12; rm2l m21; rm22 m22; rxc=xc; rdx I dxl; rdx2 =dx2; ryc yc; rdylI dylI; rdy2 dy2; for (k kv~h0;k; kkOhO rvOhOlik]; if rxc [kOhO] xmin rxc[kOhO] xmin; else if rxc[kOhO] xmax rxc[kOhO] xmnax; if ryc[kOhO] ymin ryc[kOhO] ymin; else if ryc[kOhO] ymnax ryc[kOhO] =ymax; xcr rxc[kOhO] dn; ycr ryc[kOhO] dn; mx floor((double)xcr); my floor((double)ycr); rdx I k~h0] xcr mx; rdy 1 [kOhO] ycr my; *::*rdx2 [kOhO] 1.0 -rdx [kOhO]; rdy2[kOhOI 1.0: rdyl[kOhO]; ::mx =kx +dn *mx; =ky +dn *my; *::rmllI[kOhO my *nx +mx; rml 2[kOhO] rm 1[kOhO]; *if (dxlkOhO] 0.0) rm2[kOhO] =rml2[kOhO] dn; rin2l[kOhO] rml I[kOhO]; :.if(dyl[kOhO] !=0.0)rm2l[kOhO] =rm2l[kOhO] +nxtdn; rm.22[kOhO] rml[kOhO]; if (dxl[kOhO] !=0.0)m22[kOhO] =rm22[kOhO] +dn; if dyl1[kOhO] =0.0 rm22[kOhO] rm22[kOhO] nxtdn; ej1~k (n =(dd +dd)) nxtdp nx dd dd; nxtdrn nx dd dd; x t ln=dd-kx;xmax=nx-dd-kx-1; ymin dd ky; ymax fly dd ky 1; (register long kc, kOhO, *rmll1, *lrn12, *rm.2; register float *rxc, *rdxl, *rclxa; register float *ry(c, *rdyl, *rdy2; rmlIl ml rm12 m12; rr2I m2l; rin22 m22; rxc =xc; rdx I dx1; rdx2 =dx2; ryc =yc; rdy I dy1; rdy2 =dy2; for kkv~h0; k; kkOhO if rxc[kOhiO] xmin rxc[kOhO] xmin; else if rxc[kOhOI xmax).-xc[kOhO] xmax; if ryc[kOhO] ymin ryc[kOhQ] ymnin; else if rycikOhO] yrnax ryc[kOhO] ymax; -131xcr (ryc[kOhOI rxc[kOhO]) /dn; ycr (ryc[kOhO] rxc[kOhO]) /dn; mx floor((double)xcr); my floor((double)ycr); rdxlI[kOhO] xcr -mx; rdy 1 [kOhO] ycr my; rdx2[kOhO] 1.0 rdxl[kOhO]; rdy2[kOhO] 1.0 rdy 1[kOhO]; mxx kx dd my); myy ky dd *(mx my); rm I I[kOhO myy *nx +mxx, rmlI2[kOhO] rml I kOhO]; if dx I[kOhOI 0.0 r I 2[kOhO] rm Il2[kOhOJ nxtdp; rm2l [kOhO] rm I [kOhO]; if dy 1[kOhO] 0.0 rm2l.[kOhO] rm2I [kOhO] nxtdm; rm22[kOhO] rm I I[kOhOJ; if dxlI kOh0] 0.0 rm2')2[kOhO] rm22 [kOhO] nxtdp; if (dyl[kOhOI!= 0.0 )rm22[kOhO] nn22[kOhO nxtdm; return (0) retum(1); wanI(dd, dii, mnh, my, nh) dd, dii, nih, my, nh; *o

S

registe- long i, j, ix, jy, ji1; long nhtdn, nhpl, nhml; nhtdn nh dii; nhpl nh +1; nhnil nh -1; if di dd) kv~h0 my mh kvphp (my 1) (nih 1); kv~hp =mv (nih kvmhp =(mv 1) *(mh 1); kvnihm=(mv *(mh kvCim mv 4" (nih 1); kvphm=(mv 1)*(mh 1); kvph (mv mh; )else if dn (dd dd))( kv~h0 my nih (mv kvphp (mv 1) *(nih 1) (my 1) 1); kv~hp my *C pmh +(mv (rdi-2); kvmhp (my 1) (nih 1) (mv 1) (nih 1); kvnmh0= (mv *mh (mv *(mh kv~hrnimv (mh +(mv (nh kvphm=(mv 1) (mh 1) +(mv 1) (ih- 1); kvph (mv nh (mv (nh )else refturn(1); i ;j 1; for (jy 0; jy mv; jy++) -132i =jI;jil +=nhtdn; for ix 0; ix mh; ix++) mvOhO[i] j i++;ji dn; i 1;jl =1; for (jy It;jy mv;jy++) j jl; jl nhtdn; for ix 1; ix mh; ix++) mvmhm j mvphp[i] j nhpl I dd; i4; j dn; i= 1;jl 1; for (jy 0; jy mxv; jy++) j jl; jl nhtdn; for (ix 1; ix <mh; ix++) mv~hm~i] =j nv~hp[i] j dn; j dn; i=1;jl=dn*nh+l; for (jy=1; jy <mv; jy++) j. j j1 +--nhtdn; for ix( 1; ix <m;i+ mvphxn[i] j; 'a.mvn'hp[i] j nhml *dd; j dn; i=1; j1 dn *nh 1; for (jy= 1; jy <mv; jy++) j =j1; jlI nhtdn; for ix 0; ix mh; ix-i+) mvphO[i] j mvmhO~i] j nhtdn; i++;ji dn; if dn dd.) return(0), i niv n 1; j I= dd nh +dd +1; -133for (jy 1; jy <mv; jy+-+ j =j1; j I nhtdn; for (ix ix <mh; ix-H mvOhO[i] j j dn; 1 (my 1) (mh 1) 1; i dd nh +dd +1; for (jy= 1; jy <mv; jy++) j =j1; j I +=.nhtdn; for (ix 1; ix <mh; ix++) mvmhm~i] j; mvphp[i] j nhpl *dd; dn; mv (mh +1; i I jldd nh +dd +1; for (jy 1; jy <mv; jy++) j =j 1;j I nhtdn; for (ix ix <mh; ix++) rnmv~hm[i] j mvOhp[i] j dii; j dn; jl =dd nh, o-dd +1; for (jy 1; jy <mv; jy++) j j 1; jlI ahtdn; for (ix 1; ix <mh; ix+-imvphm~i] j mvmhp~i] j nhml I dd; j dn; i=(my 1) *mrh+ 1; jl =(dn+dd)*nh+dd+ 1; for (jy jy <mv; jy++) j =jl1;j I nhtdn; for (ix 1; ix <mh; ix++) mvphO[i] j mvmhO[i] j nhtdn;, j dn; -134return(O); conjgr(rknf, feps, kmax) long kmax; float *irinf, feps; register long k, kOhO; long kit, ksv, ksc, nsc; float aiphak, betak, rezd, fnsp, t-ps; float *ru, *pu, *qu, *zij, *du, *ud, *aj; float *rv, *qv *zv, *dv, *vd, *bi; long *mll1, *mf2, *n121, *m22; float *w1 1, *w12, *w21, *w22; double ak, akml 1.0, bk, dtemp, sqrto; dtemp 0.0; for ksp 0; ksp nsp; ksp++) ru =arsd[ksp].rucv; rv arsd[ksp].rvcv; Pu =arsd[ksp].pucv; pv arsd[ksp].pvcv; qu =arsdfksp].qucv; qv arsdfksp].qvcv; du =arsd[ksp].duvf;, dv arsd[ksp].dvvf; =arsd[ksp].udvf; d=alkp.df for k kvOh0;,k; kkOhO =mv~h0[k]; rulikOhOl =du[kOhO]; pufkOhO] =udlikOhO]; rv[0I dv[Oj; i.:pv[0]=vd[0I; rnattvro; nsc arsd[ksp].cns; if nsp I1) nsc for ksc 0; ksc nsc; s+ ksv arsd[ksp].cnv[ksc]; ml.. Iu arsd[ksp].mll[ksc]; m12 arsd[ksp].nil2[ksc]; m21. arsd[ksp].mn2l [ksc]; m22 arsd[ksp].m22[ksc]; wi 1 arsd[ksp].w 1[ksc]; w12 arsd[ksp].w1I2[ksc); w21, arsd[ksp].w21[ksc]; w22 arsd[ksp].w22[ksc]; ud arsd[ksv].udvf;, vd arsdliksv].vdvf;, for (k kv~hO; k; kkOhO =mnv~h0[k]; qu[kOhO] w 11[kOhO] ud[m I I[kOhO]] wl2[kOhO] ud[ml2[kOhOI] w2l [kQ0] ud[m2l [kOhO]] w22[kOhO] ud[ni22[kOhO]J; qv[0] kv~h0 axsd[ksp].wrnn[ksc] vd[0]; -135ru =arsd[kspj.rucv; rv arsd[ksplj.rvcv; qu =arsd[ksp].qucv; qv =arsd[ksp].qvcv; zu =arsd[ksp].zucv; zv =arsd[ksp].zvcv; ai =arsd[ksplj.ainv; bi =arsdfksp].binv; for k kv~h0; k; rl< jh0I rukOhO] qu[kOhO]; ZL'-t 3] ailikOhO] ru[kOhO]; dtemp =dtemp u[kOhO] rujikOhO]; rv[0] rv[0] qv[O]; zv[O] bi[0] rv[0]; Itemp dtemp rv[0] rv[O]; *rknf sqrt(dtemp); if C*rknf feps return; fnsp nsp; teps fepsffnsp; Vfor (ksp ksp <nsp; ksp++) {uas~s]rc;r rdkp.vv pu =arsd[ksp].rucv; v arsd[ksp.vcv; Pu =arsd[ksp].pucv; pv arsd[ksp].pvcv; qu =arsd[kspl.qucv; qv arsdllksp].qvcv; arsd[ksp].zucvf; v arsd[ksp].zvcv;a: u arsd[kspj.aidvf; asd[kspl.dvf; for kit 0; kit <kmax; kit++) *dtemp=0.0; for kkv~hO; k; kOhO =mv~h0[k]; :dtemp dtemp ru[k~h0] ru[kOhO]; dtemp dternp rvilO] r[] :rczd sqrt(dtemp); if (rezd teps) break; kOhO niv~hO[k]; A A -u[k~h~j zu[kOhO]; Ak Ak rv[0] zv[OJ; if kit> 0) betak A/ alanl; for kkvh; k; kkOhO =mvh[k]; pulikOhO] pu[kOhO] betak zu[kOhO]; pv[0] O~v[0] betak zv[0]; )else([ for k kv~h0; k; k-) k mO=nv~h0[k]; -136pu[kOhO] zu[kOhO]; pv[0] =zv[O]; xnattvro; bk 0.0; for k kv~h0; k;

I

kOhO =mv~h0[k]; ba.= bk pu[kOhO] *qu~kOhO];

I

bk bk pv[0I qvIIO); aiphak =A bk; for kkvOh0; k; kkOhO; mvOh0[k]; ud[kOhO] alphak *pu~kOhO]; *ru[kOhO] -=alpbak ukh] zu[kOhOj= ai[kOhO] ru[kOhOj; vd[0] aiphak pv[0]; ry[0] -=aiphak qv[O]; zv[0] bi[0] rv[0]; aimI ak; strntx(dd, dn, nh, nv, ri) long dd, dn, nh, nv; float *rkf, *:::*legister long k, mn, lcml, km2; long s, ki, *ml, *mp2, itst; Sdouble dutnp, dvtmp, dutdu, temnp, tmnpl, tinp3, sqrtO, expo; -T float *sumap, *sinud, *uuvf, *uvf, *auul, *bvvl, scn; float *sumb, *smvd, *vwvf, *vtvf, *dvvf, *cuvl, *ajrw, *binv; Coat *gg, *aun, scssq, scs, scd; ijuvf Prsd[ksp].uuvf; vvvf arsd[ksp].vvvf, iitvf arsd[ksp].utvf; vtvf arsd[ksp].vtvf;, duvf =arsd~ksp].duvf;, dvvf arsd[ksp].dvvf;, sumna arsd[ksp].rucv; sumb =arsd[ksp].rvcv; smud arsd[kspl.zucv; smvd =hisd[ksp].zvcv; amy ars~I[ksp].ainv; binv arsd[ksp].binv; auu I arsdlikspl.auul1; bvvl1 arsd[ksp].bvvl1; cuvI= arsd[ksp].cuvl; suinb[Q] 0.0; snivd[0] 0.0; for( k =kv~h0; k; kin rnv~h0fk]; surna[in] gui; smud(in] gun; sumb[0] gvl; sm. vd[0] gvn; if ist setvtr(suma, sumb, smud, smvd, dd, dn, nh, 0) -137fprintf(stderr,"Subroutine setmtx: setvtr fprintf(stderr,"returnes value equal to %dMn', itst); return(itst); scn dn; sod sqrt(two) dd; scn 1.0 sn; scd 1.0 sd; for 9 2; s 6; switch (s) case 2: k1 =kvphp; ml I mvmhm; mf2 =mvphp; gg =arsdfksp].g2sq; auu =arsd[kspl.auu2; O scs =scd; break; case 3: k I=kv~hp; ml =mvOhm, m2 =mv~hp; gg =arsd[ksp].g3sq; *.ai =arsd[ksp].auu3; break; case 4: ki I kvuihp; ml I mvphm; m2 =mvnihp; gg =arsd[k.V].g4sq; :auu arsd[ksp].auu4; break; case k1 =kvnih0; =mvphO; m2 =mvmh0; ::gg auu scs =son; break; scssq scs sos; for (k kl;kk--) kmIn mlI kxn2 =m2I[kl; dutrnp scs (uuvf[km2] uuvf[kml]); dutdu =dutrnp *duunp; temp (arsd[ksp].csq arsd[ksp].bsq gg[km2]) dutdu; tinpI arsd[ksp].a scssq exp((double)(-1.0 temp)); tmp3 arsd[ksp].a scssq exp((double)(-3.0 temp)); auu[kxn2] tmp3; suma[kml] trnp3; suma[km2] tmp3; duvf[kmnl] uuvf[km2] tmpl;,, sxmud[kml] tmpl; duvf[km22] uuvf~lkml] tmpl; srtlud[km2I tnpl; -138temp 0.0; for( k =kvOh; k; m mvOh0[kl; auu 1 suma[m]; duvff ml utvffm] smud[m] uuvf[m]; temp temp dluvf[m] duvf[m]; sumb[0]; dvvflOl vtvf[O] smvd[0] vvvf[0I[ temp temp dvvf[0] dvvfl0l; *rkf *irkf sqrt(temp); for (k =kv~h0; k; k--

I

m mv~h0[k]; ainv[m] 1.0 /auI[] binv[0] 1.0 bvl[01; return(0); egister long k; float *pu, *pv, *qu, *qv; :..float *av~h0, *avphp, *av~hp, *avmhdp, *avmJh0, *cvOhO; &Oh0 arsd[ksp].auul; cv~h0 arsd[ksp].cuvl; wphp arsd[ksp].auu2; av~hp arsd[ksp].auu3; !;avmp arsdfksp].auu4; avmhO Pu =arsd[kspl.pucv; pv =arsd[ksp].pvcv; =a arsd[kspj.qucv; qv =arsd[kspll.qvcv; .(register long kOhO; =arsd[kspi~bvl[0] pv[0]; for kkv~h0; k; k-- ~~h0 mv~h0[kI; qu[kOhO] avOh0[kOhOl pulkOhO] cv~h0[kOhO] *pv[0]; qvO]= qv[0] cvOh0[kOhO] pu[kOhO]; iegister long kmhm, kphp; for( k =kvphp; k; kmhm mvmhim[k]; kphp mvphp[k]; qu[kmhm] qu[kmhm] avphp[kphp] pu[kphpl; (register long kOhmn, k~hp; for( k =kv~hp;kk--) k~hm mvOhni[k]; k~hp mvOhp[k]; qullk~hm] qu[kOhm] av~hp[k~hp] pu[k~hpl; (register long kphm, kmhp; for k kvmhp; k; k-) -139kphmn rnvphm[k); kmnhp mvrnhpfkl; qu[kphrn] qu[kphmil avrnhp[kinhp] pu[kmnhp); (register long kphO, kmhO; for k =kvmhO; k; k-) kphO =mvphO[k]; kinhO rnvmhO[k]; qutkphO] =qu[kphOl avmhO[krnhO] pu[kmhO]; (register long kmhm, kphp; for k kvmnhm; k; k-) kmhm mvmhmllk]; kphp mvphp[kI; qu[kphp] qu[kphp] avphp[kphp] pu[kmhm]; (register long kOhin, kOhp; for k kvOhm; k; k- 0* kOhm mvOhm[k]; kOhp vhk] qu[kOhp] qu[kOhp] avOhpllkOhpll pu[kOhm]; (register long kphmn, kmhp; for k kvph; k; k-) **,.kphm =mvphm[k]; kxhp =mvinhp[k]; ::.qu[kmhp] qu[kmhpl avrnhp[kznhp] pu[kphm]; :~register long kphO, kmhO; (k =kPhOk;k-) kphO =mrvpbO[k]; kmhO =mvmhO[k]; qu[kmhO] =qu[kmhO avmhO[kmhiO] *pu[kphO]; S&ind(hc, vc, ml 1, m12, m2l, mn22, Ahl, dvi, dh2, dv2, dd, dn, nh, nv) long *mlI 1, *ml2, *rm21, *m2 dd, dii, nh, fly; gloat *hc, *vc, *dhl, *dvl, *dh2, *dv2; This subroutine perform the following functions. 1* 1. The input vector field (hc[mvOh')"[k]],vc[mvOhO[k]I) I k 1, ,kvOhO) f* is xansformed into the current vector field 1* ((hc[mvOhO[k]],vc[mvOhO[k]]) I k 1, ,kv~h0) using the relations 1* hc[niv~hO[k]] niax(l.0,min(nhhc[mv~hO[kII))

P**

1* vc[mvOhO[k]] max(1.0,rnin(nvvc[mvOh0[k]])) 1.0; 1 14 0for k 1* 2. Integer valued fields (ml ltmvOhO[k]I, m12[mvOh0[k]I, m2l[mv~h0[kJ], m22[mv~h0[k]]), k kv~h0 are determined such that the metrix 1* ml 1[mv~hO[k]) ml2[mvOhO[k]I m22[mvOhO[k]] /~forms a grid cell of size dn containing the vector 1* (hcllmvohO[k] vc[mvOh0[k]]), for every k 1, ,kv~h0. 1* 3. Real valued fields (dxl dyl1[mv~h0[k]], dx2[mv~hO[kII, dy2[mvOh0[k]]), k 1..,kv~h0 are determined such that the function value 1* f(xc[mv~h0[k]],yc[mv~h0[k]]) is ai bilinear interpolation of /~its values at the grid points 1* f(ml 1[mv~h0[k]]) f(ml2[mv~h0[k]]) f(m2l1[mvOh0[k]]) f(m22[mvOh0[k]]) defined by the relations L:f(xc[mvOh0[k]], yc[mv~h0[k]]) 1* dy2[mvOh0[k]] (dx2[mv~hO[k]] f(ml. I dx I[mv~hO[k]] f(m 12[mv~h0[k]]) 1* dIl[mv~hO[k]l dx2[mv~h0[k]] f(m2l dxl[mv~hO[k]] f(m22[mv~h0[k]])). float fnh, fnv; register long k, ktmp; double flooro, rths, r-tvs, fdn; long trnpl, tmp3, rthm, rtvm, nhmlI; long trnp2, tmp4, rthp, rtvp, nvml; nhml nh 1; nvml nv- 1; fnh nh Il; fnv nvmlI; fdn dn; for( k =kvOh; kk--) ktmnp mv~h0[k]; hc[ktmp] if hc[ktrnp] 0.0 hc[ktnip] 0.0; if hc[ktmp] fnh hc [ktmp] fnh; vcllktmp] if vc[ktmpl 0.0 vc[ktmp] 0.0;

M

-141if vc [ktrnp] fnv vc [ktmp] fnv; if Oni dd) for( k =kvOhO; k; knnp mv~h0[k]; rths hc[ktmp] fdn; rthrn floor(rths); rt's vc[ktmp] fdn; rtvm floor(rtvs); dh I[ktmp] rths rthm; dh2[ktrnp] 1.0 dhlI ktrnp]; dv 1[ktmp] rtvs rtvm; dv2[ktmp] 1.0 dv I[ktmpl; rthm rthni dn; rtvm rtvm n rthp =rthm; if (rthp <nhml rthp =rthp +dn; rtvp =rtvm; if (rtvp <nvml rrvp =rtvp +dn; m 11l[ktmp] nh *rtvm rthm 1; m12[ktmp nh rtvm +rthp 1; m2l [ktmp] nh rtvp +rthm 1; m22[ktmp] nlirtvp +rthp 1; if((rthm <0)II(rthrn nhml)U(rthp O)II(rthp nhml)) fprintf(stderr, ktrnp, hcllktmp]); fprintf(stderr, "dh ",ktmp, dhl1[ktmpD); fprintf(stderr, "dh2[%d] f, ktmp, dh2[ktrnp]); fprintf(stderr, t 'rthm rthm); fprintf(stderr, t "rthp rthp); return(1), *if((rtvxn 0)II(rtvm nvrnl)II(rtvp <0O)II(rtvp nvml)) fprintf(stderr, ktmp, vc[ktmp]); fprintf(stderr, "dvl[%d] ktmnp, dvl[ktmp]); fprintf(stderr, 'dv2[%d] ktmp, dv2[ktmp]); fprintf(stderr, "rrvm rtvm); fprintf(stderr, "rvp rtvp); return(1); return(0); It (dn (dd dd)) for( k =kv~h0; k; kktmp mvOh0[k]; rths (vc[ktmnp] hc[,ktmp]) /fdn; rtvs (vc~ktmp] hc[ktmp]) fdn; rthm floor(rths); rtvrn floor(rtvs); dh I[ktrnpl rths rthm; dh2[ktrnp] 1.0 dhlI ktrnp]; dvlI[knmp] rtvs -rtvm; dv2[ktmnp] 1.0 -dv 1[ktmp]; rthp rthm rtvm; nvp rthrn rrvm; rthm rhp frti= rp dd; tnpI= nh *rtvm rthm 1; tmp2 =rh (rtvm +dcn) +rthm 1; trnp3 =nh (rtvm +dd) +rthm +dd 1; tmp4 =nh (rtvm +dd) +rthm -dd +1; -142if (rtvm -dd ml 1 [krrnp tmp2; else m 11[ktmp] tipl1; if (rtvm nvm 1) ml2[ktmp] rrpl; else if (rthm nhml) ml2[ktmp] =tmp4; erse m12[ktmpl t).p3; if (rtvm ==nvmlI) ni2l [ktmp] =tmpl1; else if rthm 0 m21 [ktrnp] =tmnp3; else m2l [ktmp] mp4; if (rvm>= nvml dd) m22[ktmp] tmnpl; else 19W m22[ktmp] tmp2; if (rthn 0) 11 (rthm nhm1))

I

fprintf(stderr, "ch: dd fprintf(stderr, "nhml ",nhml); .fprintf(stdlerr, "rthm rthm); return(l); *fprintf(stderr, "cv: dd dd); fprintf(stderr, "nvml nvmi); return(l); return(0); return(1); setvtr(suma, sumb, smud, smvd, dd, dn, nh, nv) float *sunla, *sumb, *smud, *smvd; long dd, dii, nh, nv;

I

register long 1, k, ktnip; double dtrnp, expO; long itst, nsc, ksc, ksv, *ml 1, *m12, *nrpl, *mp22; float tmp". ruOO, tull1, tul2, tu21, tu22, twl1, tw2l, twl2, tw22, one float *lltvt, *uuvf, *hc, *dhl1, *dh2, *w 11, *w 12, *wOO, temp, dit, wxy, nhf;, float *vtvtf, *vwF, *vc, *dl *dv2, *w21, *w22, *dusq, wmn, wuv, nvf; urvf arsd[ksplutvf;, hc wrk05; nhf nh vtvf arsd(ksp].vtvf-; vc wrk06; nvf nv for( k =kv~h0; k; k--

I

ktrnp mv~h0[k]; utvflktmp] gun arsdfksp].u~vftktnipl; -143vtvfTO] kvOhO gvn arsd[kspbvOvf[O]; nsc arsdfksp.I.cns; if nsp <2 11 (nsc 1 )retun(O); for ksc 0; ksc nsc; ksc++)

I

wmn =arsd[ksp].wmn[ksc]; wxy =arsd[ksp].wxy[ksc] scm sgmw; wuV arsd[ksp].wuv[ksc] scm sgmw; dlt arsd[ksp].dltiksc]; ksv arsd[ksp].cnv[kscl; d Ii arsdfksp].wllI[ksc]; dh2 arsd[ksp].w 12[kscl; dvi arsd[ksp].w2l~ksc]; dv2 arsd[ksp].w22[ksc]; wi 1 arsd[ksp].wl 1 [ksc]; w 12 arsd[ksplj.wl2[ksc]; w21 arsd[ksp].w2l[ksc]; w22 arsd[kspl.w22[ksc]; ml 1 arsd[ksp].mll[ksc]; m12 arsd[ksp].ml2[ksc]; m2l arsd[ksp].m21 [ksc]; m22 arsd[ksp].m22[ksc]; uuvf arsd[kspl.uuvf; vvvf arsd[ksp].vvvf; for( k =kvOh0; k; ktmp vhk] hc[ktmp] hv~hO[ktrnp) clt *uuvf[ktmp); vc[krrnp] vv~hOllktmp] dlt vvvf[0]; if((itst setind(hc,vcml 1 ,ml2,m2l1,m22,dhl1,dvl1,dh2,dv2,dd,dn,nh,nv)) ***,fprintf(stdlerr, "Subroutine setvtr. setind :fprintf(stdlerr,"retuxnes value equal to itst); return(itst); dusq arsd[ksv].dusq; for kkvOh0; k; :ktnip mv~h0[k]; tuOO arsd[ksp].uuvf~ktnipj; tull =arsd[ksv].uuvffmlljktmp]]; tmpt ucl tuQO; dtmp wxy wuv dusq[ml11[ktnip]]) trnpu trnpu; wi 1 I wmn dv2[ktmp] dh2[ktmp] exp(dtmp); tu12 =arsd[ksv].uuvflml2[ktmnp]]; trnpu =tu12 u; dtznp wxy wuv dusq[m12[ktp]]) mpu tmpu; tw12 =wrnn dv2[ktrnp] dhlI[ktmp] exp(dtmp); tu2l =arsd[ksv].uuvf[m2l [ktmnp]]; tmpu =tu2l tuQO; dtmp wxy wuv dusq[m21[ktrnp]]) tpu tmpu; tw21 wnin dvl[ktmnpl dh2[ktmp] exp(dtnip); tu22 =arsd[ksv].uuvf~m22[ktmpj]; tmpu =tu22 tuQO; dtmp wxy wuv dusq[m22[ktmp]]) *tmpu tmnpu; tw22 =wmn dvlI[ktmp] dhlI[ktmp] exp(dtmp); if (hvOh0[ktnip] one) 11 (hvOh0[ktmp] nhf) 11 (vvOh0[ktinp] one) 11 (vv~hO[ktmpl nvf)) twl 1 wmn dv2[ktmpl dh2[ktmp]; tw12 wmn dv2[ktmp] dhl[ktmp]; tw2l wmn dvlllktmnp] dh2[ktmp]; tw22 wmn. dvlI[ktrpD] Ad1 [ktrnp]; temp 0.0; tmpu 0.0; wl I[ktmp] twl1; temp +=twll1; tmpu tullI twl1; -144w 12 [ktrnp] twl12; temp twl12; tmpu tul12 tw 12; w2 I[ ktrnp] tw2 1; temp tw2l1; tmpu tu2l tw2l1; w22[ktrnp]j tw22; temp tw22; tmpu tu22 tw22; suma[ktmp] temnp; smud[ktmp] temp; urvflktrnp] trnpu; suznb[O] kvOhO *wmn; sm'%7d[O] kvOhO *wmnr vtvffO] kvOhO *wmn *arsd[ksv].vvvflOl; return(O);

Claims

1. A method, comprising the steps of: capturir-q pair of stereo images of a scene separated in a horizontal direction; and creating an intermediate image between the pair by estimation of velocity vector fields constrained in the horizontal direction.

2. A method as recited in claim 1, further comprising producing a depth image from the pair and the intermediate image.

3. A method as recited in claim 1 or claim 2, wherein step creates plural intermediate images.

4. A method as recited in any one of claims 1 to 3, wherein step includes aligning an entirety of the pair and the intermediate image in a vertical direction simultaneously. A method as recited in any one of claims 1 to 4, wherein step (b) comprises the steps of: determining correspondence between the pair by constraining correspondence searching in the horizontal direction; and (ii) interpolating between the pair using the correspondences. An apparatus for producing a depth image, comprising: 20 a capture device capturing a pair of stereo images; a computer connected to said capture device and creating an intermediate image by horizontally constrained correspondence interpolation and creating the depth image from the pair and the intermediate image; and a depth image output device connected to said computer and outputting 25 the depth image.

7. A method according to claim 1 substantially as herein described with reference to the drawings.

8. An apparatus for producing a depth image substantially as herein described with reference to the drawings. DATED: 25 March, 1997 PHILLIPS ORMONDE FITZPATRICK RA, Attorneys For: ^g EASTMAN KODAK COMPANY V A ibC 'VINWORANDREASPIEClaI 6Cn.DO C re I I 66,645 ABSTRACT OF THE DISCLOSURE The present iivention performs an interpolation operation between a pair of stereo images to produce plural intermediate images. The interpiolation operation involves estimating the velocity vector field of the intermediate images. The estimation for a particular intermediate image includes constraining the search for the correspondences between the two images in the horizontal direction allowing the system to arrive at a more optimal solution. Excessive gap filling by smoothing is not necessary. At the same time the process vertically aligns the entirety of the images. .e e I