Devices for refractive field visualization

Abstract

An apparatus according to an embodiment of the present invention enables measurement and visualization of a refractive field such as a fluid. An embodiment device obtains video captured by a video camera with an imaging plane. Representations of apparent motions in the video are correlated to determine actual motions of the refractive field. A textured background of the scene can be modeled as stationary, with a refractive field translating between background and video camera. This approach offers multiple advantages over conventional fluid flow visualization, including an ability to use ordinary video equipment outside a laboratory without particle injection. Even natural backgrounds can be used, and fluid motion can be distinguished from refraction changes. Embodiments can render refractive flow visualizations for augmented reality, wearable devices, and video microscopes.

Claims

1. An augmented reality apparatus comprising: a video camera with an imaging plane, the video camera configured to capture video of a scene including a textured background, the video including captured light that passes from the textured background to the imaging plane through a refractive field; and an augmented reality display configured to render the video of the scene with a visualization of actual motions of the refractive field, the actual motions determined by correlating, over time, representations of apparent motions of the textured background observed at the imaging plane in the video from frame to frame by modelling translations of the refractive field as causing the apparent motions.

2. The apparatus of claim 1, wherein the visualization includes representations of velocities of the refractive field based upon the correlated representations of the apparent motions.

3. The apparatus of claim 1, wherein the refractive field is a fluid.

4. The apparatus of claim 1, further including a processor configured to determine the actual motions by correlating the representations of the apparent motions through an assumption of frame-to-frame changes in at least one of intensity and phase resulting from the motions alone.

5. The apparatus of claim 1, further including a processor configured to determine the actual motions by correlating the representations of the apparent motions using an approximation that the apparent motions are constant from one video frame to another video frame.

6. The apparatus of claim 1, wherein the visualization of actual motions of the refractive field is a visualization of a fluid turbulence, updraft, downdraft, or vortex.

7. The apparatus of claim 1, wherein the actual motions of the refractive field arise from motion of one or more refractive objects in a fluid, and wherein the representations of the apparent motions are representations of motions of the one or more refractive objects.

8. The apparatus of claim 7, wherein the one or more refractive objects include at least one of a smoke particle, molecule, region of the fluid differing in temperature from a surrounding region, region of fluid differing in pressure from a surrounding region, droplet, bubble, exhaust, hydrocarbon, or volatile substance.

9. The apparatus of claim 1, wherein the textured background is a natural background.

10. The apparatus of claim 1, wherein the display is wearable and further configured to render the video of the scene with a visualization of actual motions of the refractive field for viewing by a wearer's eye.

11. The apparatus of claim 10, further comprising a processor configured to determine the actual motions.

12. The apparatus of claim 10, wherein the video camera and the augmented reality display form part of an air traffic control system, an onboard airplane system, a weather monitoring system, a leak detection system, a civil engineering system, an aeronautical engineering system, or a ballistics or combustion monitoring system.

13. A wearable device comprising: a video camera with an imaging plane, the video camera configured to capture video of a scene including a textured background, the video including captured light that passes from the textured background to the imaging plane through a refractive field; and an augmented reality display configured to render a visualization of actual motions of the refractive field in a wearer's eye, the actual motions determined by correlating, over time, representations of apparent motions of the textured background observed at the imaging plane in the video from frame to frame by modelling translations of the refractive field as causing the apparent motions.

14. A video microscope comprising: a video camera with a microscopic field of view and an imaging plane, the video camera configured to output video by capturing light that passes from a textured background to the imaging plane through a refractive field in the microscopic field of view; and a display configured to show a visualization of actual motions of the refractive field in the microscopic field of view, the actual motions determined by correlating, over time, representations of apparent motions observed at the imaging plane in the video from frame to frame by modelling translations of the refractive field as causing the apparent motions.

15. The video microscope of claim 14, wherein the visualization includes representations of velocities of the refractive field based upon the correlated representations of the apparent motions.

16. The video microscope of claim 14, wherein the refractive field is a fluid.

17. The video microscope of claim 14, further including a processor configured to determine the actual motions by correlating the representations of the apparent motions using an assumption of frame-to-frame changes in at least one of intensity and phase resulting from the motions alone.

18. The video microscope of claim 14, further including a processor configured to determine the actual motions by correlating the representations of the apparent motions using an approximation that the apparent motions are constant from one video frame to another video frame.

19. The video microscope of claim 14, wherein the motions of the refractive field arise from motion of one or more refractive objects in a fluid, and wherein the representations of the apparent motions are representations of motions of the one or more refractive objects.

20. The apparatus of claim 19, wherein the one or more refractive objects include at least one of a molecule, region of the fluid differing in temperature from a surrounding region, region of fluid differing in pressure from a surrounding region, droplet, bubble, hydrocarbon, or volatile substance.

21. The video microscope of claim 14, wherein the textured background is a natural background.

22. The video microscope of claim 14, wherein the textured background is an artificial textured background having microscopic features.

23. The video microscope of claim 14, wherein the microscopic field of view includes at least a partial view of a chemical sample or a cell culture or other biological sample.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

(2) The foregoing will be apparent from the following more particular description of example embodiments of the invention, as illustrated in the accompanying drawings in which like reference characters refer to the same parts throughout the different views. The drawings are not necessarily to scale, emphasis instead being placed upon illustrating embodiments of the present invention.

(3) FIG. 1A is an image-based overview of refractive flow principles, including a comparison of results based on optical flow and those based on refractive flow.

(4) FIG. 1B is an overview of refractive flow similar to the overview in FIG. 1A, but with different target scenes and textured backgrounds.

(5) FIG. 1C(a)-(f) illustrates use of refractive flow according to embodiment methods and apparatus to obtain velocity and depth of a refractive field.

(6) FIG. 1D(a)-(d) illustrates results of using uncertainty measurements to improve refractive flow measurements.

(7) FIG. 2 illustrates use of an embodiment of the invention for an augmented reality display for an airplane pilot.

(8) FIG. 3A illustrates use of an embodiment device in air traffic control.

(9) FIG. 3B illustrates use of an embodiment in a wearable device, namely augmented reality eyewear.

(10) FIGS. 3C-3E illustrate various environments in which embodiments can be used for monitoring applications, including wind farms, civil engineering, and hydrocarbon leak monitoring.

(11) FIG. 3F illustrates an embodiment apparatus that can be used to monitor a chemical process.

(12) FIG. 3G illustrates and embodiment apparatus including a video microscope to monitor a biological cell culture.

(13) FIG. 4A is a block diagram illustrating an embodiment device.

(14) FIG. 4B is a block diagram illustrating an embodiment apparatus that includes a display and an alarm signal.

(15) FIG. 4C is a block diagram illustrating an embodiment apparatus for determining fluid flow in three dimensions.

(16) FIGS. 5A and 5B illustrate example procedures for correlating representations of motions in a video and using the correlated representations.

(17) FIG. 5C illustrates an example procedure for determining uncertainty in correlated representations of motion.

(18) FIG. 6A illustrates use of a brightness constancy principle in the prior art to track motion of a solid object.

(19) FIG. 6B illustrates use of a refractive constancy principle according to embodiments of the invention to obtain fluid object flow or refractive flow.

(20) FIG. 6C illustrates use of a refractive constancy principle according to embodiments of the invention to obtain depth of a refractive field as well as refractive flow.

(21) FIG. 7 illustrates the results of refractive flow analysis using a single camera to view various scenes.

(22) FIGS. 8A-8E illustrate quantitative evaluation of a refractive flow method according to an embodiment of the invention.

(23) FIGS. 9A-9D illustrate a controlled experiment performed to evaluate fluid flow velocities estimated using refractive flow method according to an embodiment of the invention.

(24) FIG. 10 is illustrates results of fluid depth measurements based on stereo video sequences.

(25) FIGS. 11A-11B illustrate the effect of background texturing on depth measurement results with and without uncertainty applied.

(26) FIG. 12 illustrates quantities referred to in proving refractive constancy.

(27) FIG. 13 illustrates quantities used in proving stereo refractive constancy.

DETAILED DESCRIPTION

(28) A description of example embodiments of the invention follows.

(29) The teachings of all patents, published applications and references cited herein are incorporated by reference in their entirety.

(30) Measuring and visualizing how air and fluid move is of great importance, and has been the subject of study in broad areas of science and technology, including aeronautical engineering, combustion research, and ballistics. Special techniques may be necessary to visualize fluid flow, and an example is shown in FIG. 1A.

(31) FIG. 1A shows images 102 of a burning candle 104. A plume 105a of the candle 104 is invisible to the naked eye and to standard video techniques. The heat rising from the burning candle 104 causes small, invisible distortions of a textured background 106 seen behind the candle 104. The small distortions of the background result from light rays being refracted by the plume 105a as they travel from the background 106 to a camera (not shown). As the air in the candle plume 105a moves, small changes in the refractive properties of the air appeared as small visual deformations, or motions, of the textured background 106. These motions can be referred to as wiggles or refraction wiggles.

(32) Multiple techniques have been proposed for visualizing or measuring such small distortions in fluids, such as sound tomography, Doppler LIDAR and schlieren photography. Current techniques to visualize and measure fluid flow can be divided into two categories: those which introduce a foreign substance (dye, smoke, or particles) into the fluid; and those that use the optical refractive index of the fluid.

(33) Where foreign substances or markers are introduced, the motions of the markers can be used to give a visualization of the fluid flow. Through particle image velocimetry (PIV), a quantitative measurement of the flow can be recovered by tracking the particles introduced into the fluid.

(34) In the optical technique category, schlieren photography is one method. Schlieren photography uses a calibrated, elaborate optical setup that can include a source grid placed behind a photography subject, a cutoff grid incorporated into a special camera, and special lamps for illumination to amplify deflections of light rays due to refraction to make the fluid flow visible. To attempt to remedy some of the limitations of schlieren photography, the schlieren method has been simplified in the past decade. In background-oriented schlieren (BOS) photography (also referred to as synthetic schlieren), the complicated, hard-to-calibrate optical setup of traditional schlieren imaging is replaced by a combination of a video of a fluid flow captured in front of an in-focus textured background (such as the textured background 106 shown in FIG. 1A) and optical flow calculations (further described hereinafter). The refraction due to the fluid motion is recovered by computing optical flow between each frame of the video and an undistorted reference frame. A resulting displacement field gives a visualization of the fluid flow.

(35) Continuing to refer to FIG. 1A, this figure provides an example of schlieren photography images, namely the center images 110. The synthetic schlieren photography images 110 were created by applying optical flow calculations to the raw input video frames 102. The application of optical flow calculations is symbolized in FIG. 1A by the arrow 103a. The resulting displacement field provides a visualization of the fluid flow above the candle 104. The displacement field can be presented visually in various ways, including the visible plume 105b in the schlieren photography images 110 and the apparent motion velocity vectors v(x,y) 114 (or representations of the motions or refraction wiggles) overlaid on the video frames 112. The velocity vectors v(x,y) 114 are oriented in an xy plane with an x-axis 117 and y-axis 119.

(36) However, previous techniques have multiple disadvantages. Introduction of foreign substances or markers with use of PIV requires additional expense and is not suitable for continuous monitoring or virtual reality applications, for example. Even for optical techniques, setups can be expensive and limited to laboratory use. A schlieren photography setup, as described above, is an example of such limitations.

(37) BOS presents an additional issue, in that a textured background must be in focus. To address this issue, light-field BOS has been proposed. In light-field BOS, a textured background like the textured background 106 can be replaced by a light field probe in certain instances. However, the light field probe comes at a cost of having to build a light field probe as large as the fluid flow of interest. Besides, for many applications such as continuous monitoring and virtual reality, light field probes would be undesirable or even impossible to use in some cases.

(38) Furthermore, simply computing optical flow in a video will not generally yield the correct motions of refractive elements, because the assumptions involved in these computations do not hold for refractive fluids. Optical techniques such as schlieren photography rely on a brightness constancy assumption. Under the assumption of brightness constancy (explained further hereinafter), any changes in pixel intensities are assumed to be caused only by translation of the photographed objects. However, even if the brightness constancy holds for solid, reflecting objects, brightness constancy does not typically hold for refractive fluids such as the plume 105a-b of the candle 104 in FIG. 1A. In the case of refractive fluids, changes in pixel intensities result from both translation and refraction of the fluid. Thus, optical flow does not yield actual motion of a refractive fluid, but rather how much an image projected through the refractive fluid onto the imaging plane of a camera moves.

(39) Disclosed herein is a novel method to measure and visualize fluid flow. Even complicated airflows, for example, can be measured, without elaborate optical equipment. Binocular videos are not necessarily required, and two-dimensional (2D) fluid flow can be derived even from standard monocular videos. Embodiments disclosed herein apply a recognition that local changes in intensity may be locally invariant, and an observed scene can be modeled as being composed of a static background occluded by refractive elements that bend light and move through space. Embodiments of the current invention include computing the optical flow of the optical flow by re-deriving the optical flow equation for a proposed image formation model. As used herein, video denotes any format that includes a series of images of a scene over time, and video camera denotes any instrument capable of capturing such a series of images. Video cameras can be camera configured to capture a series of images at a standard 30 frames per second (fps), or a higher or lower speeds, or a camera that takes single photographs at similar rates, not necessarily at a constant frame rate.

(40) Embodiments can rely on fluid flows of interest containing moving refractive structures, which are almost always present in fluid flows of interest. Embodiment methods can track features that move with the fluid to recover the fluid flow (also called refractive flow, or calculated or actual fluid flow). Advantageously, fluid flow can be obtained without the added cost of seeding a fluid with particles and are thus markerless. Embodiments can use the displacement field (the field resulting from an initial computation of optical flow) to construct features to track and can be implemented outside a laboratory and even on mobile devices.

(41) Continuing with the description of FIG. 1A, the velocity vectors v(x,y) 114 are referred to as apparent motion velocity vectors because, as described above, simply applying optical flow based on brightness constancy fails due to detected intensities being caused by both translation of the plume 105 and refraction imparted by the plume. Changes in the apparent motion velocity vectors between the frames 112a, b, c . . . can be used as an indicator of the motion of the transparent fluid, and the disparity of these motion features between viewpoints of two cameras, for example, can be used to determine the depth of the fluid. Optical flow may be applied again, as represented by an arrow 103b, to obtain refractive flow according to embodiments of the invention. Optical flow is applied again by a processor, which is configured to correlate, over time, the representations of motions in the video frames 112 from frame to frame as a function of motion of the refractive field of the plume 105.

(42) Correlating representations of motion can be done in accordance with methods described hereinafter in conjunction with FIGS. 5A and 5B, for example. Correlating over time can include tracking or matching observed, apparent motions represented by, for example, the vectors 114, over the video frames 112a, b, c, and so forth, for example. Representations of motions can additionally be correlated over space by tracking the representations over video frames from a single camera or multiple cameras. The correlation can be done as a function of motion of a refractive field where, in the correlation, the represented, observed are assumed to arise from motion of a refractive field, with a background of the refractive field being stationary, for example.

(43) Light from the textured background 106 passes through the candle plume 105 to an imaging plane of a video camera (not shown). Optical flow is applied to obtain apparent motion velocity vectors v(x,y) 114, and the apparent motion velocity vectors of the video frames 112 are correlated through applying optical flow again to obtain refractive flow.

(44) The results of refractive flow are shown in FIG. 1A in the form of calculated fluid flow velocity vectors u(x,y) 118 overlaid on video frames 116. Together, the video frames 116 with the calculated fluid flow velocity vectors u(x,y) 118 constitute a representation of a scene, which varies over time in successive video frames 116a, 116b, 116c, and so forth.

(45) FIG. 1B is a refractive flow overview 120, similar to the refractive flow overview 100 in FIG. 1A, showing results of applying example methods described hereinafter in fluid flow applications. The applications illustrated in FIG. 1B include wind visualization 125a, landing 125b of a plane 124, and takeoff 125c of the plane 124. Each application includes a different textured background 126. In the wind visualization 125a, the textured background 126a consists primarily of a natural background including trees. In contrast, the textured backgrounds 126b and 126c contain many elements made by humans. Thus, FIG. 1B illustrates how a wide variety of different textured backgrounds can be used, including both natural and artificial elements, to obtain refractive flow measurements with standard video images.

(46) In FIG. 1B, the process of obtaining refractive flow is similar to that illustrated in FIG. 1A. Namely, optical flow is applied to raw input video frames 122 to obtain representations of apparent motion observed at the imaging plane of the camera (not shown). These representations are presented in the form of video frames 132 with apparent motion vectors v(x,y) 134 overlaid the representations in the video frames 132 are correlated by again applying optical flow to obtain refractive flow measurements, which are represented by the video frames 136 with calculated fluid flow vectors u(x,y) 138 overlaid.

(47) FIGS. 1C(a)-(f) illustrate application of refractive flow principles to obtain both fluid flow and depth of a refractive field. FIG. 1C(a) shows a video image obtained from a left-side video camera, and FIG. 1C(b) shows the same scene as in FIG. 1C(a), but from a right-side camera (not shown). Heat rising from two burning candles 104 causes small distortions of the background due to light rays refracting as they travel from the background to the cameras, passing through the hot air created by the candles. FIGS. 1C(c) and (d) show visualizations of the heat disturbances obtained by the existing synthetic schlieren imaging for the left- and right-side cameras, respectively. The synthetic schlieren images in FIGS. 1C(c) and (d) can be used to visualize the small disturbances and reveal the heat plumes, but schlieren imaging is unable to produce a measurement of the actual motion of the heat plumes. Under reasonable conditions described hereinafter, the refraction patterns (or observed motions) move coherently with a refracting fluid, allowing one to measure accurately both two-dimensional motion of the fluid flow and depth of the fluid flow.

(48) FIG. 1C(e) shows calculated fluid flow velocity vectors 118 that can be obtained from a monocular video. In this case, the vectors 118 were obtained from images from the left-side camera. Further, as shown in FIG. 1C(f), the depth of the fluid flow can be obtained from a stereo image sequence. FIG. 1C(f) shows a disparity between the position of the plume obtained from the left-side camera and the position obtained from the right-side camera, and the disparity can be used to determine the depth of the refractive field.

(49) FIGS. 1D(a)-(d) illustrate how refractive flow measurements can be further improved by applying uncertainty principles described more fully hereinafter. FIG. 1D(a) shows a raw image of the candles 104, and FIG. 1D(b1) shows apparent motion velocity vectors obtained from applying optical flow. The apparent motion of a fluid that is visualized using optical flow can also be described as a wiggle. FIG. 1D(b2) shows a schlieren image produced according to existing methods. FIG. 1D(c) illustrates apparent motion velocity vectors depicting the heat plumes rising from the candles 104, and the apparent motion velocity vectors are obtained using methods and devices according to embodiments of the invention. FIG. 1D(d) shows the apparent motion velocity vectors weighted, as described hereinafter more fully, to improve accuracy refractive flow results even further.

(50) FIG. 2 illustrates an application of embodiments of the invention that includes real-time visualization of turbulence that can be encountered by an airplane 224. Disturbances such as the updraft 228a and the turbulent airflow to 228b create a refractive field 228 that could be encountered by the airplane 224. Light rays from a textured background 226 pass through, and are refracted by, the updraft 228a and turbulence 228b and are captured by a video camera 240 mounted on the airplane. Also within the field of view 242 of the video camera are various elements of the textured background 226, including a landscape 226a, horizon 226b, cloud 226c, and distant turbulent air 226d. Other elements of the textured background 226 can be useful at night, including city lights 226e, moon 226f, and stars 226g. Thus, FIG. 2 illustrates that a wide variety of different background elements can contribute to a textured background.

(51) Even where turbulent air 228b in the refractive field 228 is being monitored and visualized in the foreground of the camera 240, for example, distant background turbulent air 226d is relatively constant in position compared with the foreground refractive field 228 and is useful as part of the textured background 226. The video camera 240 provides raw video images 244 to be stored in memory 246. A fluid flow processor 248 obtains the raw video images 244 and applies example methods described hereinafter to produce refractive flow measurements. The refractive flow measurements are output from the fluid flow processor 248 in the form of images 228 with calculated refractive flow vectors overlaid. The images 228 are presented in an display 250 to be seen by a pilot 252. The display 250 is updated continually to provide the pilot 252 with a real-time view of the refractive field 228, and thus the pilot 252 can see and attempt to avoid turbulence such as the updraft 228a and the turbulent air 228b.

(52) It should be noted that in other embodiments, the display 250 can be an augmented reality display. For example, calculated fluid flow velocity vectors such as the vectors 118 in FIG. 1A can be projected onto the windows of the airplane 224 so that the vectors are seen in conjunction with the visual view of the atmosphere in which the refractive field 228 is located.

(53) FIG. 3A illustrates another setting in which embodiments of the invention can be used, namely in an air traffic control context. A video camera 340a is mounted on an air traffic control tower 354 and has a video camera field of view 342a. The field of view 342a includes areas in which airplanes, such as an airplane 324, take off and land. In other embodiments, the video camera 340a is mounted near runways or in other locations where it is useful to monitor airflow.

(54) As understood in aviation, airplanes can produce wake turbulence, such as jetwash 328a and wingtip vortex 328b produced by the airplane 324. Images from the video camera 340a can be used to provide a display 350a for an air traffic controller 352. Moreover, the video camera 340 can see other air disturbances, such as wind 329. The embodiment of FIG. 3A thus illustrates how safety can be improved by using embodiments of the invention to monitor in real-time for both conditions created by airplanes and also natural weather conditions that can affect aviation. A real-time display, such as the display 350a, can be a much more useful to tool to identify these conditions than windsocks or the other weather monitoring instruments currently used because the display 350a can show a wide field of view and is not limited to point locations as windsocks are.

(55) FIG. 3B illustrates how embodiments of the invention can be used in conjunction with augmented reality eyewear 356. A video camera 340b is mounted to the eyewear 356. The video camera 340b has a field of view 342b. An augmented reality display image 350b is generated in the wearer's eye by a projector (not shown). The display image 350b includes an image of the field of view 342b with a representation of fluid flow in the field overlaid on the image. The representation of the fluid flow can include calculated fluid flow velocity vectors, such as the vectors 138 shown in FIG. 1B. However, in other embodiments, fluid flow is represented in an augmented reality display image by color coding or some other means, for example.

(56) Augmented reality eyewear, such as the eyewear 356 shown in FIG. 3B, can be useful in many applications in which it is useful to visualize fluid flow. Such situations can include hang gliding, for example, wherein a real-time knowledge of updrafts, downdrafts, and other air currents is important for safety and flight optimization. It should be noted that alternative embodiments of the invention can include virtual reality eyewear in which a virtual reality display image is generated in the wearer's eye. In the case of virtual reality, the eyewear would be opaque, and a completely artificial image would be generated in the wearer's eye. This is in contrast to the augmented reality eyewear 356 in FIG. 3B, in which the wearer can see outside the eyewear, and the view is simply augmented by an overlaid image 350b. According to embodiments of the invention, either fixed video cameras such as the camera 340a in FIG. 3A or mobile cameras with eyewear like the eyewear 356 in FIG. 3B can be used to visualize airflow or refractive flow in a variety of situations. Some of these situations are shown in FIGS. 3C-3E.

(57) FIG. 3C illustrates a wind farm 358 in which wind turbines 360 are placed. Fixed or mobile video cameras can be used to obtain video images from the wind farm 358 to optimize output of wind turbines 360. For example, the wind turbines 360 can be positioned to maximize output based on refractive flow around each turbine.

(58) FIG. 3D illustrates a city 362 in which there are various refractive elements that can be visualized by embodiment apparatuses. Such refractive sources include a steam source 364a, an exhaust source 364b, a heat and particulate matter source 364c, and wind 364d from a water body 366 near the city 362. A knowledge of such refractive sources can be used in a variety of ways in civil engineering and environmental monitoring applications, for example.

(59) FIG. 4C is a block diagram illustrating a correlating apparatus 400c that can receive stereo representations 474d of motions captured by two video cameras and process the representations to obtain a depth 430 of the refractive field 428. The light 427 from the textured background 426 passes through the moving refractive field 428 and on to the imaging plane 470 of the video camera 440 and an additional imaging plane 470 of an additional video camera 440 private. The video camera 440 outputs representations 440c of motions in the captured video, and the video camera 440 outputs representations 474c of the motions in the video from camera 440. Both the representations 474c and 474c are included in the stereo representations 474d of the motions in the videos. The stereo representations 470d are stored in a correlating apparatus 400c in a memory 446c. The correlating apparatus 400c also includes a 3-D (or depth) processor 448c, which obtains the stereo representations 474d from the memory 446c. The 3-D processor 448c is configured to correlate, over space, the stereo representations 474d in the two videos from frame to frame as a function of the motion 429a of the refractive field 428. Correlation of the stereo representations over space can include matching an apparent motion feature observed I one location of the imaging plane 470 with the same apparent motion feature observed in another location of the imaging plan 470. The 3-D processor 448c also calculates, based on the correlated representations of the motions, a depth 430 of the refractive field. The 3-D processor 448c can optionally output depths for multiple points on the refractive field 428.

(60) The 3-D processor 448c can optionally output calculated 2-D fluid flow velocity vectors 438. An optional uncertainty processor 448d calculates uncertainties of fluid flow velocity measurements and depth measurements according to embodiment methods and weights the measurements to output weighted fluid flow velocity measurements and depth measurements 432. While the uncertainty processor 448d is separate from the 3-D processor 448c in FIG. 4C, in other embodiments, a single processor performs correlation, depth, and uncertainty calculations. The uncertainty processor 448d can be configured to determine an uncertainty in the representations of the motion or in the depth of the refractive field using a variety of techniques described hereinafter. In particular, determining the uncertainty can include calculating a variance of the representations of motion and weighting the representations of motion as a function of the variance. Furthermore, the processor 448d can weight the representations of motion as a function of a degree of texturing of the textured background 426. In particular, in regions of greater texturing, the measurements of motion, or representations of motion, can be given greater weight. The processor 448d can apply an L2 norm to the representation of motion or calculate a covariance of the concatenation of a plurality of representations of motion, as described hereinafter. The determination of uncertainty can include waiting the representations 474c and 474c as a function of a logarithm of a covariance of the representations of motion, and the covariance can be covariance of the representations 474c from video camera 440 and the additional representations 474c from the video camera 440. Furthermore, in other embodiments that include display such as the display 450 in FIG. 4B, the representations of motion can be displayed with the uncertainties in the representations by any means available, including, for example, shading and color coding the representations of motion.

(61) FIG. 3E illustrates one example application for the petroleum industry. Namely, stationary or mobile cameras, such as those shown in FIGS. 3A-3B, can be used to monitor a pipeline 366 for a hydrocarbon leak 368. The pipeline 366 can be in a refinery or in a drilling location or other field or distribution location, for example. Furthermore, any volatile substance, droplet, bubble, or other refractive substance or object can be visualized where it differs in refractive properties from a surrounding fluid.

(62) FIGS. 3F and 3G illustrate that embodiments of the invention are not limited to outdoor or even to macroscopic refractive flow applications. FIG. 3F shows an embodiment of the invention used for monitoring a chemical process 351. The chemical process 351 occurs in a chemical processor structure 349 and is monitored by a video camera 340b configured to acquire video images of the chemical process 351 occurring in a fluid. Images of the chemical process 351 are sent to a fluid flow processor 248, which correlates representations of fluid motions in the chemical process 351 to obtain refractive flow measurements. The fluid flow processor 248 provides measurements of fluid flow in the chemical processor 349 and sends the measurements to a display 350c for viewing by a worker 352.

(63) Refractive flow in the chemical processor 349 can result from bubbles rising, fluid heating, chemical reactions, or any other process that produces fluid currents within the processor 349. The video camera 340b is configured to capture light passing from an artificial textured background 326 through a fluid in the chemical processor 349 to an imaging plane (not shown) of the video camera 340b. In a controlled, artificial environment such as the chemical processor 349 in FIG. 3F, the artificial textured background 326 is preferable because it is provides superior texturing for visualizing refractive changes in the chemical processor 349 and is much easier to implement in a controlled setting than in the outdoors. However, in other embodiments, a natural textured background, such as the textured background 126a in FIG. 1B, can be used.

(64) FIG. 3G illustrates an embodiment of the invention used to view fluid motion in a biological cell culture 353. FIG. 3G illustrates that embodiments of the invention can be used to visualize fluid flow even on a microscopic scale. The biological cell culture 353 is situated in the video microscope 343 and is viewed by a microscope objective lens 341, which is part of the video microscope 343. The microscope 343 also includes a video camera 340c that acquires images of the biological cell culture 353 viewed through the microscope objective lens 341 and captured at an imaging plane (not shown). A textured background (not shown) can include, for example, an artificial textured background such as the background 326 in FIG. 3F. The textured background can be scaled appropriately to have features comparable in size to the fluid motion in the biological cell culture 353. The video camera 340c outputs images of the cell culture 353, or representations of motions observed in the cell culture 353, to the fluid flow processor 248. The processor 248 correlates the representations as a function of motion of the refractive field in the cell culture 353. The correlated representations of the motions in the fluid are provided to a display 350c for viewing as calculated fluid flow vectors or in any other convenient form.

(65) Lighting (not shown) in the video microscope 343 can be adjusted to optimize illumination of the artificial textured background. As understood in the art of microscopy, either forward lighting or backlighting can be used, as deemed appropriate. It should also be noted that, although not shown in FIGS. 3F and 3G, dual video cameras could be used in conjunction with the devices shown to obtain depth measurements of the refractive field in the chemical processor 349 or in the biological cell culture 353, as previously described in conjunction with FIGS. 1C(a)-(f), and as more fully described hereinafter.

(66) FIG. 4A illustrates a correlating apparatus 400a including memory 446a and a correlation processor 448a. The memory 446a is configured to store representations 474a of motions in a video. The video is obtained from a video camera 440 with an imaging plane 470. The correlating apparatus 400a also includes a processor 448a that is configured to correlate, over time, the representations 474a of the motions in the video that are stored in the memory 446a. The processor 448a is configured to correlate the representations from frame to frame as a function of motion 429a of a refractive field 428 through which light 427 from a textured background 426 passes to the imaging plane 470 of the video camera 440. One way to perform this correlation is to apply the two-step refractive flow method described hereinafter. The correlation processor 448a is configured to output correlated representations 418. The correlated representations 418 are an optional example output, and the correlation processor 448a outputs other types of outputs in other embodiments. The correlated representations 418 can be output in the form of calculated fluid flow velocity vectors, such as the velocity vectors 118 in FIG. 1A, for example.

(67) In FIG. 4A, the representations 474a are in the form of electronically stored apparent motion velocity vectors v(x,y) stored for various times for different video frames, such as the vectors 114 in FIG. 1A. However, in the other embodiments, the representations stored in the memory can be in any form that permits the processor to correlate, over time, the representations from frame to frame as a function of motion of the refractive field. These representations can include raw video images, for example. Either of those forms, and other forms of representations of the motions not shown here, can be suitable for correlation of the representations of the motions by the processor 448a with appropriate adjustments to correlation methods. In particular, a variety of electronic formats can be used to store the images showing the motions as understood in digital photography. Furthermore, as described earlier in conjunction with FIG. 1A, and as more fully described hereinafter, one convenient form of representations of motions in the video includes apparent motion velocity vectors such as the vectors 114 shown in FIG. 1A and obtained by initially applying optical flow. Moreover, while the memory 446a is shown as part of the correlating apparatus 400a in FIG. 4A, in other embodiments, the memory 446a is incorporated into the video camera 440.

(68) The correlation processor 448a is configured to output correlated representations 418 of motions in the video. The correlated representations 418 are an example output from the processor and can be in the form of calculated fluid flow velocity vectors such as the vectors 118 in FIG. 1A, for example.

(69) FIG. 4B is a block diagram illustrating an alternative correlating apparatus 400b. The apparatus 400b includes several optional features not included in the apparatus 400a in FIG. 4A. The memory 446b and processor for 448b in correlating apparatus 400b also are configured to operate somewhat differently from those in the apparatus 400a of FIG. 4A. The memory 446b in FIG. 4B for 446b and processor 448b are part of a computer 472. The memory 446b is configured to receive raw video images 444 from the video camera 440. Raw video frames 402 are transferred from the memory 446b to the processor 448b for performing initial optical flow calculations. The processor 448b then transfers representations 474b of motions in the raw video to the memory 446b in the form of apparent motion velocity vectors, such as those represented graphically in FIG. 1A as vectors 114. The processor 448b then obtains the apparent motion velocity vectors 114 from the memory 446b to apply optical flow calculations again to obtain refractive flow measurements. Refractive flow measurements are then output by the processor 448b to a display 450 in the form of calculated fluid flow velocity vectors 438 overlaid on an image of the scene to present an augmented reality scene 416.

(70) The output from the processor 448b to the display 450 can be via a dedicated graphics card, for example. The processor 448b also outputs an alarm signal 476 to an alarm 478. The alarm signal 476 can be output to the alarm 478 via a driver (not shown), for example. The alarm 478 is activated by the alarm signal 476 when at least one calculated fluid flow velocity vector exceeds a given threshold in magnitude. Thus, as the alarm 478 feature demonstrates, displays such as the display 450 having augmented reality scenes are not the only way to use refractive flow measurements. The refractive flow measurements can be, optionally, further processed, displayed, or used in any way desired to visualize or calculate or represent fluid flow.

(71) FIG. 4B also illustrates that the refractive field 428 generally includes many motions 429b, 429c, and 429d, for example, throughout the field. The motions 429b-d of the refractive field can arise from motion of one or more refractive objects in a fluid. Furthermore, the representations 474b of motions in the raw video can be representations of motions of the one or more refractive objects. Refractive objects can include smoke particles, molecules, regions of fluid differing in temperature from a surrounding region, regions of a fluid differing in pressure from surrounding region, droplets such as rain or condensation, bubbles such as bubbles in a liquid, exhaust, hydrocarbons, volatile substances, or any other refractive elements or substances in a fluid.

(72) The 3-D processor 448c is configured to correlate, over space, the representations 474c and 474c of motions from the two respective video cameras 440 and 440 with the respective imaging planes 470 and 470. The 3-D processor 448c correlates the representations over space from frame to frame as a function of the motion 429a of the refractive field 428. The representations of motions 474c and 474c are correlated over space when the location of a representation 474c observed at the imaging plane 470 is matched to a location on the imaging planes 470 where a corresponding representation 474c is observed. Based on the correlated representations of the motions, a depth of the refractive field 428, or a point thereon, can be calculated, as described more fully hereinafter.

(73) FIG. 5A illustrates a procedure 500a for obtaining fluid flow or refractive flow measurements. At 580a, video captured by a video camera with an imaging plane is obtained. At 580b, representations of motions in the video are correlated over time by a processor. The representations are correlated from frame to frame as a function of motion of a refractive field through which light from a textured background passes to the imaging plane of the video camera.

(74) FIG. 5B illustrates an alternative procedure 500b for obtaining fluid flow measurements. Referring to FIG. 5B, at 582a, video captured by a video camera with an imaging plane is obtained. At 582b, representations of motions in the video obtained are defined to be representations of apparent motions observed at the imaging plane. The representations of the apparent motions are also determined by assuming that frame-to-frame changes in at least one of the intensity and phase are caused by the motions alone. At 582c, the representations of motions in the video are correlated over time from frame to frame as a function of motion of a refractive field through which light from a textured background passes.

(75) In FIG. 5B, the correlation at 582c includes sub-elements 582c1, 582c2, and 582c3. At 582c1, the textured background is modeled as a stationary background, and the refractive field is modeled as translating between the textured background and the video camera. At 582c, the apparent motions are considered to be constant from one video frame to another video frame. At 582c3, actual motions of the refractive field are estimated based on the determined representations of the apparent motions by assuming that motions of the refractive field are an exclusive cause of the apparent motions.

(76) Referring to FIG. 5B, at 582d, a velocity of the refractive field is calculated based upon the correlated representations of the motions. At 582e, a fluid velocity vector field is calculated, based on the correlating, as a function of position in the refractive field.

(77) At 582g, a representation of the fluid velocity vector field is displayed. The display can be like the display 450 in FIG. 4B, for example. At 582h, the representation of the fluid velocity vector field is displayed in an arrangement with the real-life scene to present an augmented reality scene. The augmented reality scene can be in the form of video frames like the video frames 116 with the calculated fluid flow velocity vectors 118, as shown in FIG. 1A, for example. Alternatively, the augmented reality scene can be in the form of images such as the images 228 shown in FIG. 2 on the display 250 or such as the augmented reality display image 350b generated in the wearer's eye in FIG. 3B. Furthermore, in other embodiments, a virtual reality scene can be generated instead of an augmented reality scene.

(78) Continuing to refer to FIG. 5B, at 584, the fluid velocity vector field is used to monitor for a hydrocarbon leak, a fluid turbulence, updraft, downdraft, or vortex. As shown in FIG. 5B, 584 can occur after calculating the fluid velocity vector field at 582e. However, as also shown by a dashed line, in other embodiments, the fluid velocity vector field is used for monitoring at 584 after displaying an augmented reality scene at 582h. This difference in when the fluid velocity vector field is used for monitoring relates to whether or not an augmented reality scene is used for display purposes. As explained in conjunction with FIG. 4B, the fluid velocity vector field may be used in an augmented reality display such as the display 450, or the fluid velocity vector field may be used only for further calculations or monitoring performed by a computer, for example.

(79) One example of monitoring performed by a computer is found in FIG. 4B, wherein the processor 448b analyzes the fluid velocity vector field to determine whether any velocity exceeds a given threshold. If the threshold is exceeded, the processor sends an alarm signal. In this way, use of the fluid velocity vector field can be by means of a computer and/or further calculations and need not be used to present an augmented reality or virtual reality scene for viewing by a person.

(80) FIG. 5C illustrates a procedure 500c for determining an uncertainty in refractive flow measurements. At 581a, a representation of motion in a video captured by video camera is modeled as arising from motion of a refractive field through which light passes from a stationary background to an imaging plane of the video camera. At 581b, uncertainty in the representation of the motion is determined by a processor.

(81) Determining the uncertainty in the representation can include calculating variance of the representation of motion and weighting the representation of motion as a function of the variance. In some embodiments, determining the uncertainty includes weighting the representation of motion as a function of the degree of texturing of the stationary background. For example, where the background is heavily textured, the wiggles or representations of motion observed by a video camera can be determined with greater certainty. On the other hand, where a background is only lightly textured, the representations of motion can be determined with less certainty, and these representations can be downweighted in determining refractive flow. Furthermore, where the representations of motion or refractive flow measurements determined there from are displayed on a monitor or other image, uncertainty in the representations can also be displayed through, for example, shading or color coding the representations of motion. As described more fully hereinafter, determining the uncertainty can include applying an L2 norm to the representation of motion, calculating a covariance of the concatenation of a plurality of representations of motion, or weighting the representation of motion as a function of the logarithm of a covariance of the representation of motion in the video captured by the video camera and an additional representation of the motion in an additional video captured by an additional video camera, such as the video cameras 440 and 440 in FIG. 4C.

(82) Refractive Fluid Flow

(83) An example goal of refractive flow is to recover the projected 2D velocity, u(x, y, t), of a refractive fluid such as hot air or gas, from an observed image sequence, I(x, y, t). Before proceeding by analogy to a differential analysis of refraction flow, a differential analysis of standard optical flow for solid objects is first presented hereinafter.

(84) Differential Analysis of Optical Flow

(85) FIG. 6A includes a diagram 600a illustrating application of a brightness constancy assumption to flow of a solid object 627 according to known methods. As seen in FIG. 6A, the solid object 627 moves with time from one position at time t.sub.1 to another position at time t.sub.2. An object point 629a on the solid object 627 also moves between times t.sub.1 and t.sub.2. The object point 629a is imaged to a projection point 629a on the camera imaging plane 670 at times t.sub.1 and t.sub.2, with a center of projection 690 for the video camera. An apparent motion velocity vector 614a represents the motion of the object point 629a between times t.sub.1 and t.sub.2 as observed at the camera imaging plane 670. Because light rays do not penetrate the solid object 627, the textured background 426, shown for purposes of comparison with other figures, is not relevant in FIG. 6A.

(86) Under the brightness constancy assumption, any changes in pixel intensities, I, are assumed to be caused by a translation with velocity v=(v.sub.x, v.sub.y) over spatial horizontal or vertical positions, x and y, of the image intensities, where v.sub.x and v.sub.y are the x and y components of the velocity, respectively. That is,
I(x,y,t)=I(x+x,y+y,t+t),(1)
where (x, y) denotes the displacement of the point (x, y) at time t+t.

(87) Applying a first order Taylor approximation to Equation 1 gives

(88) $\begin{array}{cc}\frac{I}{x}{v}_x+\frac{I}{y}{v}_y+\frac{I}{t}=0.& \left(2\right)\end{array}$

(89) The velocity v can then be solved for by using optical flow methods such as the Lucas-Kanade or Horn-schunck methods known in the art.

(90) Extension to Refraction Flow

(91) Brightness constancy, however, does not hold for refractive fluids, as rays passing through a point on the fluid intersect the background at different points at different times, as shown in FIG. 6B.

(92) FIG. 6B includes a diagram 600b illustrating application of a refractive constancy principle to flow of a refractive fluid 628 according to embodiments of the invention. In contrast to a solid object 627, the refractive fluid 628 can transmit at least some light from the textured background 426 through the refractive fluid 628 to the camera imaging plane 670. Thus, as described previously, motion observed at the camera imaging plane 670 arises both from motion of the refractive fluid 628 and from refraction of light at the refractive fluid 628 during its motion. The textured background 426 can be a natural background, for example, in that the textured background 426 can include natural or everyday objects and need not be a special laboratory setup. For example, the textured background 426 can be like the textured backgrounds 126A-C in FIG. 1B or light the textured background 226 in FIG. 2, for example.

(93) The refractive fluid 628 is shown at two time instants t.sub.1 and t.sub.2. At time t.sub.1, the fluid 628 refracts light rays, effectively moving background points 631 and 633 on the background to different positions, or projection points 629b, on the camera imaging plane 670. At time t.sub.1, for example, the background point 631 appears at a different position on the imaging plane 670 than it does it at the time t.sub.1+t. Similarly, at time t.sub.2, the background point 633 appears at a different point on the camera imaging plane 670 than it does at time t.sub.2+t.

(94) In FIG. 6B, the observed motions on the imaging plane are referred to as apparent motion velocity vectors 614b1 and 614b2. In general, the direction of the observed motions on the imaging plane due to the refraction can be arbitrary. In the case illustrated in FIG. 6B, the observed motions or apparent motions 614b1 and 614b2 are in the same direction as the actual or calculated fluid flow, as shown by a calculated fluid flow velocity vector 638. However, the direction of the observed motions on the imaging plane due to refraction can be consistent over time as the fluid changes its position. According to embodiments of the invention, the observed motions can be tracked over time to recover a projected 2-D motion of the refractive fluid 628.

(95) To obtain a calculated fluid flow such as the calculated fluid flow velocity vector 638 shown in FIG. 6B, the observed scene can be modeled as composed of a stationary background, (x, y) (such as the textured background 426), and a refractive field, r(x, y, t), (such as the refractive fluid 628) translating in between the scene and the camera:
I(x+r.sub.x(x,y,t),y+r.sub.y(x,y,t)={tilde over (I)}(x,y),(3)

(96) The refractive field, r, indicates the spatial offset in the imaging plane, caused by light rays bending due to index of refraction gradients between the scene and the camera. For example, such refractive gradients can be due to bubbles of hot air traveling with the velocity of moving air.

(97) The inventors have recognized that it can be further assumed that the moving refraction gradients imaged to a point (x, y) have a single velocity, u(x, y). In other words, the apparent motions v(x,y) are considered to be constant from one video frame to another frame. Under these conditions, namely when motions, or translation, of the refractive field is assumed to be the exclusive cause of local image velocity changes, this refraction constancy can be exploited to explain the observed local velocity changes by the translation, u(x, y), of an assumed refractive field. In other words, the frame-to-frame motions in a video are considered to be a function of motion of the refractive field r(x,y,t). This yields a refraction constancy equation:
r(x,y,t)=r(x+u.sub.xt,y+u.sub.yt,t+t).(4)

(98) It is shown hereinafter in the section labeled Proof of v=r/t that, under the condition of Equation 4, the temporal derivatives of the refraction field correspond directly to the observed motion in the sequence. That is, v=r/t. Additionally taking the partial derivative of Equation 4 with respect to t, the observed motion can be written as
v(x,y,t)=v(x+u.sub.xt,y+u.sub.yt,t+t).(5)

(99) Since, the v(x,y,t) constitute frame-to-frame representations of motions in a video, and since r(x,y,t) in Equation 4 is expressed as a function of motion u(x,y) of the refractive field, an implementation of Equation 5 by a processor is, thus, one way to correlate, over time, representations of motions in the video from frame to frame as a function of motion of the refractive field.

(100) It should be pointed out that the correlating performed in embodiments of the present invention should be distinguished from the correlating performed in signal processing. Further, while one way of correlating representations of motions in the video as a function of motion of a refractive field is described above, the embodiments can include other ways of correlating. For example, another way to correlate representations of motions in the video is to first filter the representations of motion by complex oriented filters, and then extract the phase change of filters response between frames.

(101) Assuming v varies smoothly and that t is small, a first order Taylor expansion can again be used to approximate Equation 5 as

(102) $\begin{array}{cc}\frac{v}{x}{u}_x+\frac{v}{y}{u}_y+\frac{v}{t}=0.& \left(6\right)\end{array}$

(103) Equation 6 can be used to develop and implement a two-step approach for solving for the fluid motion, u. The observed or apparent motions, v, can first be calculated using the standard optical flow equation (Equation 2). Then, optical flow can be applied again, this time on the estimated motions, to get the actual motion of the refractive fluid. Solving for u(x,y) for at least one point (x,y), therefore, constitutes calculating a velocity of the refractive field based upon the correlated representations v(x,y,t) of the observed motions. Further, solving for several points can be carried out to produce a fluid velocity vector field. Since the observed motions correspond directly to changes in the refraction field, and following the refraction constancy assumption, computing the optical flow of the optical flow will yield actual or refractive fluid flow.

(104) Refractive Flow Versus Optical Flow

(105) Applying standard optical flow to an observed video sequence will not yield the correct motions of the refractive fluid. To illustrate this, the refraction constancy equation (Equation 4) may be considered again. Approximating Equation 4 by the first order using a Taylor expansion yields

(106) $\begin{array}{cc}\frac{r}{x}{v}_x+\frac{r}{y}{v}_y+\frac{r}{t}=0.& \left(7\right)\end{array}$

(107) Note that the optical flow result is equal to the change of refraction field, v=r/t. Then from Equation 7, it can be observed that:

(108) $\begin{array}{cc}v=\mathrm{Du},\mathrm{where}D=\left[\begin{array}{cc}\frac{r_x}{x}& \frac{r_x}{y}\\ \frac{r_y}{x}& \frac{r_y}{y}\end{array}\right],u=\left[\begin{array}{c}{u}_x\\ u_y\end{array}\right],& \left(8\right)\end{array}$
where the matrix D in Equation 8 can be called the refraction Jacobian matrix (or the refraction matrix).

(109) Equation 8 shows that the observed motions, v, are in fact the product of the refractive matrix D and the fluid flow v. This means that, under the image formation model, the observed motions are caused by both refraction gradients and fluid motion. For example, the same observed motion can result from large refraction changes and small fluid motion, or large fluid motion and small refraction changes.

(110) The form of Equation 8 provides a few insights on refractive flow. First, it can be noted that refractive flow can be obtained where flow is not completely stable laminar flow, in which case the refraction matrix, D, would be zero. This is analogous to the aperture problem known for optical flow. Second, it can be noted that optical flow (brightness constancy) alone will not be able to disambiguate fluid motion from refraction changes. In fact, with the refraction constancy assumption (Equation 4), these two types of motions are implicitly disambiguated by attributing lower-spatial-frequency motion to fluid motion, and higher-spacial-frequency motion to refraction gradients.

(111) Refractive Flow based on Phase Analysis

(112) It is useful to note that the approach to fluid flow described above provides no restriction on how to actually compute optical flow in a refractive flow methods. For example, good results can be obtained when using local phase changes for initially characterizing the observed or apparent motions, instead of using optical flow on the intensities. It should be noted that motions due to refraction can be small, even on the order of one pixel per frame or less.

(113) Phase changes can be obtained by building a two-orientation, two-level complex steerable pyramid and applying DC-balanced temporal high pass filter to the phases (e.g., convolving with [1, 1]). This allows one to avoid having to use an undistorted reference frame, as done in conventional background-oriented schlieren techniques. This decomposition provides both amplitude and phase, which, respectively, capture the texture strength and motion.

(114) A temporal Gaussian blur can be applied to an input video sequence, and an amplitude-weighted spatial blur can be applied to the phases. These variations improve the signal-to-noise ratio of the signal significantly, allowing extraction of very subtle signals due to refraction. To deal with small camera motions, the amplitude weighted mean of the phase over a single orientation, scale, and frame can be subtracted from a phase signal.

(115) After extracting and processing the phases, feature vectors can be created by combining the phase from 19 frames in time (the frame of interest itself, 9 frames before, and 9 frames after the frame of interest). Using a wider temporal extent to compute features can provide a richer representation and can improve results relative to only tracking features between single frames.

(116) Notably, even ordinary consumer cameras can be used to acquire video images in embodiments of the invention. For example, the candle, wind, and landing sequences shown in FIGS. 1A and 1B were taken with ordinary consumer DSLR cameras. However, higher performance cameras may also be used in some applications for improved results. For example, a high-speed camera was used to film the takeoff video images shown in FIG. 1B. Specifically, the takeoff sequence was filmed at 500 frames per second.

(117) Where a background is sufficiently textured, even small differences in air temperatures can be visualized with devices implementing the example methods described above. For example, in a separate test video sequence (not shown), the flow of heated air around a person was visualized. Since the temperature difference between a person and the ambient air can be very small, the amount of refraction produced by heated air around person can likewise be small, resulting in deflections as small as 100th of a pixel. In this test case, the signal-to-noise ratio was increased by having the person running up and down 20 flights of stairs and then acquiring the video of the person in a 4 C. refrigerator using a well-textured background consisting of a multiscale wavelet noise pattern.

(118) However, embodiments of the invention can also visualize scenes without an artificial background like the multiscale wavelet noise pattern previously mentioned. In the wind, take off, and landing videos, for example, atmospheric wind patterns are visualized using the forest and buildings as the textured background to visualize changes in index of refraction at a given position in the scene.

(119) It should be noted that measuring wind speed is an important factor in weather forecasting. Embodiments of the invention allow for denser sampling of wind speed than the point measurements given by an anemometer. Thus, embodiments of the invention can improve the accuracy and range of forecasts.

(120) Proof of v=r/t

(121) If the background is considered to be stationary, then the brightness constancy assumption can be applied to fluid objects:
I(x+r.sub.x(x,y,t),y+r.sub.y(x,y,t),t)=I(x+r.sub.x(x,y,t+t),y+r.sub.y(x,y,t+t),t+t).(9)

(122) A first order approximation to the right hand side of Equation 9 yields

(123) $I\left(x+r_x\left(x,y,t\right),y+r_y\left(x,y,t\right),t\right)=I\left(x+r_x\left(x,y,t\right),y+r_y\left(x,y,t\right),t+t\right)+\frac{I}{x}\frac{r_x}{t}t+\frac{I}{x}\frac{r_y}{t}t,$

(124) which yields

(125) $\begin{array}{cc}\frac{I}{x}\frac{r_x}{t}+\frac{I}{y}\frac{r_y}{t}+\frac{I}{t}=0.& \left(10\right)\end{array}$
Therefore, the observed image motion v is equal to the temporal gradient of the refractive field r/t.
Conditions for 2D Refraction Constancy

(126) In this section, it is shown when the 2D refraction constancy assumption holds. In general, a fluid element is described by a spatiotemporally-varying index of refraction n(x, y, z, t). The 2D refraction field is determined by the gradient of the index of refraction. It can be assumed that index of refraction does not change when a fluid element moves (3D refraction constancy). That is,
n(x,y,z,t)=n(x+x,y+y,z+z,t+t).

(127) Suppose a ray traveling in this field is defined as f(s)=(f.sub.x(s), f.sub.y(s), f.sub.z(s)). When the index of refraction is close to 1 (as in air), the following equation holds

(128) $\begin{array}{cc}\frac{\mathrm{dd}}{\mathrm{ds}}=n,& \left(11\right)\end{array}$
where

(129) $d=\frac{d\stackrel{->}{r}}{\mathrm{ds}}/.Math.\frac{d\stackrel{->}{r}}{\mathrm{ds}}.Math.$
is a unit vector describing the local ray direction.

(130) If the distance between the air flow and background is much larger than the thickness of the air, then the refractive field is approximately equal to

(131) $\begin{array}{cc}r\left(x,y,t\right)=-\mathrm{hP}{}_{f\left(s\right)\text{:}{C}_{x,y,t}}n\left(f\left(s\right),t\right)\mathrm{ds}& \left(12\right)\end{array}$
where h is the distance between the air and the background, P is a projection matrix that projects the motion in 3D to 2D, and C.sub.x, y, t is the ray that hits the point (x, y) on the image plane at time t. It can be assumed that all the points on the ray path C.sub.x, y, t share the same motion.

(132) When the motion of air is small, the ray hits (x, y) on the image plane at time t passes through the similar region as the ray hits (x+v.sub.xt, y+v.sub.yt) on the image plane at time t+t. Therefore,

(133) 0${}_{f\left(s\right)\text{:}{C}_{x,y,t}}n\left(f\left(s\right),t\right)\mathrm{ds}={}_{f\left(s\right)\text{:}{C}_{x+v_{x}t,y+v_{y}t,t+t}}n\left(f\left(s\right),t+t\right)\mathrm{ds},$
Furthermore, if it is further assumed that over a short time, the motion along z-dimension is negligible compared with the distance between the air and the background, then h in Equation 12 can be considered to be constant. Therefore, r(x, y, t)=r(x+v.sub.xt, y+v.sub.yt, t+t). This proves the 2D refraction constancy.

(134) Calculations have been derived for tracking the movement of refractive fluids in a single video and for recovering a 3-D position of points on a fluid surface from stereo video sequences. Both of these calculations are based on the refractive constancy assumption previously described, namely that intensity variations over time (the wiggles or apparent motions) can be explained by the motion of a constant refraction field.

(135) Calculations for Three Dimensional (3D) Measurements and Uncertainty

(136) FIG. 6C includes a diagram 600c illustrating measurement of a depth 630D from stereo video images in accordance with embodiments of the invention. Similar to FIG. 6B, apparent motion velocity vectors at a camera imaging plane are viewed as the refractive fluid moves from a time t to time t+t plus Delta T such that apparent motion velocity vectors, or wiggles, can be observed in two different imaging planes of two different cameras and can be matched or correlated over space to recover the depth of the fluid.

(137) FIG. 6C is similar to FIG. 6B in that apparent motion velocity vectors, or wiggles, are viewed using light rays passing from a textured background 426 through a refractive fluid 628 that is in motion. However, FIG. 6C includes two cameras, a left camera with a left camera imaging plane 670.sub.L and a right camera with the right camera imaging plane 670.sub.R. Left camera has a left center of projection 690.sub.L and observes wiggles 614.sub.L. The right camera has a right center of projection 690.sub.R and sees wiggles 614.sub.R as the refractive fluid 628 moves during the small time t.

(138) The wiggle features in an input video are computed using optical flow, and then using those features to estimate the motion and depth of the fluid, by matching them across frames and viewpoints.

(139) The distortion of the background caused by the refraction is typically very small (in the order of 0.1 pixel) and therefore often hard to detect. The motion features have to be extracted carefully to overcome inherent noise in the video, and to properly deal with regions in the background that are not sufficiently textured regions, in which the extracted motions are less reliable. To address these issues, probabilistic refractive flow and stereo methods that maintain estimates of the uncertainty in the optical flow, the refractive flow, and the fluid depth are presented hereinafter.

(140) The proposed methods have several advantages over existing methods: (1) a simple apparatus that can be used outdoors or indoors, (2) they can be used to visualize and measure air flow and 3D location directly from regular videos, and (3) they maintain estimates of uncertainty.

(141) Gradients in the refractive index of air can show up as minute motions in videos. Theorems are stated hereinafter that establish the relation between those observed motions in one or more cameras and the motion and depth of refractive objects in a visual scene.

(142) In the monocular case previously shown in FIG. 6B, a video camera observes a static background 426 through a single refractive fluid layer 628. When the fluid is moving, the same point on the background may trace back to different points on the image plane due to refraction changes (see the dashed and solid lines in FIG. 6B. Such observed motions can be called refraction wiggles, or wiggles for short. When a fluid object with a non-uniform surface moves a small amount from time t.sub.1 (solid shape) to time t.sub.1+t (dashed shape), the light ray passing through a point B on the background changes due to the refraction (from solid line to dashed line), producing a visual shift on the image plane (short arrow), which is the wiggle and can be represented by the apparent motion velocity vectors 614b.sub.1 and 614b.sub.2. At time t.sub.2, the fluid moves to another location, and the projection of the movement of the fluid onto the image plane (long black arrow) can be recovered. The direction of the wiggle and the direction of the refractive flow can be different, as the wiggle direction also depends on the shape of the refractive element.

(143) Because the input in FIG. 6B is only a monocular video sequence for one camera imaging plane 670, the method can at most recover a 2D projection of the motion of the refractive fluid 628.

(144) The task of refractive flow is to recover the projected 2D fluid motion from an input video sequence.

(145) Similarly, in the stereo case, two cameras with image planes 690.sub.L and 690.sub.R image static background 426 through a refractive fluid layer 628. If a standard stereo matching directly on the input stereo sequence is applied, the method will recover the depth of (solid) points in the background. In contrast, in refractive stereo it is desirable to recover the depth of the (transparent) fluid layer, stereo fusing on the motion wiggles rather than the image intensities.

(146) Methods herein after presented for both refractive flow and refractive stereo are based on this observation: for a point on the fluid object, its refraction wiggle is constant over a short time and across different views. To state this observation more formally, first the wiggle that a point on a fluid layer generates can be defined.

(147) Definition 1 (Refraction wiggle): Let A be a point on the refractive fluid layer, B be the intersection between the background and the light ray passing through A and the center of projection at time t, and t be a short time interval. Then the wiggle of A from time t to t+t is the shift of the projection of B on the image plane during this time.

(148) Then this definition leads to the following two refractive constancy theorems.

(149) Theorem 1 (Refractive flow constancy): Suppose the fluid object does not change its shape and index of refraction during a short during [t1, t2]. Then for any point on the fluid object, its wiggle v(t1) at t1 equals to its wiggle v(t2) at t2.

(150) Theorem 2 (Refractive stereo constancy): Suppose there are at least two (n2) cameras imaging a refractive fluid object, and they are all parallel and close to each other. Then at any time t, and for any point on the fluid object, the corresponding wiggle in all the cameras is the same.

(151) Proofs of Theorems 1 and 2 are included hereinafter in the Supplementary Material section.

(152) FIGS. 6B and 6C illustrate these two theorems, respectively. In the monocular case shown in FIG. 6B, when the fluid moves from time t.sub.1 to t.sub.1+t.sub.1 and from time t.sub.2 to t.sub.2+t.sub.2, it generates wiggles shown by the respective arrows, and according to Theorem 1, these two wiggles are equal. Similarly, in the stereo case in FIG. 6C, when the fluid moves from t to t+t, it generates a wiggle on the left and right views (shown by the respective arrows 614C.sub.L and 614C.sub.R), and according to Theorem 2, these two wiggles are also equal.

(153) The practical implication of these theorems is that matching the projection of a point on the fluid object across frames (over time) or viewpoints (over space) can be done by matching the observed wiggle features. That is, while the position of intensity texture features is unrelated to the fluid 3D structures, the position of wiggle features respect features on the fluid surface and can serve as an input feature for optical flow and stereo matching methods. Hereinafter, derivations of optical flow and stereo methods for tracking and localizing refractive fluids are presented.

(154) The goal of fluid motion estimation is to recover the projected 2D velocity, u(x, y, t), of a refractive fluid object from an observed image sequence, I(x, y, t). As described in the previous section, the wiggle features v(x, y, t), not the image intensities I(x, y, t), move with the refractive object. Thus, estimating the fluid's motion can be done using two steps: 1) computing the wiggle features v(x, y, t) from an input image sequence I(x, y, t), and 2) estimating the fluid motion u(x, y, t) from the wiggle features v(x, y, t). Each of these steps is described in turn below.

(155) Computing Wiggle Features

(156) The brightness constancy assumption in optical flow is that any changes in pixel intensities, I, are assumed to be caused by a translation with motion v=(v.sub.x, v.sub.y) over spatial horizontal or vertical positions, x and y, of the image intensities, where v.sub.x and v.sub.y are the x and y components of the velocity, respectively. That is,
I(x,y,t+dt)=I(x+v.sub.xdt,y+v.sub.ydt,t).

(157) Based on this brightness constancy equation, a traditional way to calculate the motion vector v is to minimize the following optical flow equation:

(158) $\tilde{v}=\arg \underset{u}{\min }\underset{x}{.Math.}{}_1{.Math.\frac{I}{x}{v}_x+\frac{I}{y}{v}_y+\frac{I}{t}.Math.}^2+{}_2\left({.Math.\frac{v}{x}.Math.}^2+{.Math.\frac{v}{y}.Math.}^2\right)$
where .sub.1 and .sub.2 are weights for the data and smoothness terms, respectively.
Estimating Fluid Motion

(159) Let u.sub.x and u.sub.y be the x and y components of the fluid's velocity. The refractive constancy equation for a single view sequence is then:
v(x,y,t+t)=v(x+u.sub.xt,y+u.sub.yt,t)

(160) Notice that refractive constancy has the exact same form as brightness constancy (Equation 1), except that the features are the wiggles v rather than the image intensities I. Applying optical flow to the wiggle features v (rather than to the intensities, I), will yield the fluid motion u. The fluid motion u by minimizing the (fluid flow) following equation:

(161) $\tilde{u}=\arg \underset{u}{\min }\underset{x}{.Math.}{}_1{.Math.\frac{\tilde{v}}{x}{u}_x+\frac{\tilde{v}}{y}{u}_y+\frac{\tilde{v}}{t}.Math.}^2+{}_2\left({.Math.\frac{u}{x}.Math.}^2+{.Math.\frac{u}{y}.Math.}^2\right)$

(162) This is similar to the Horn-Shark optical flow formulation and similar to optical flow methods, and a multi-scale iterative method can be used to solve it.

(163) FIGS. 1D(a)-(d) demonstrate a result of this refractive flow method, when applied to a video of a burning candle (candle). First, wiggle features (b.sub.1) and (b.sub.2) are extracted from the input video (a). Since wiggle features move coherently with the air flow, the method correctly estimates the motion of the thermal plume rising from the candle. Notice that the observed motions (wiggles) have arbitrary directions, yet the estimated refractive flow (c) is much more coherent.

(164) Such processing is very sensitive to noise, however, as can be seen in FIG. 1D(b.sub.1) and (b.sub.2). The problem can be even more severe for less textured backgrounds. This motivates the probabilistic formulation below.

(165) Both the refractive flow, and its uncertainty can be estimated. Consider a background that is smooth in the x direction and textured in they direction. Due to the aperture problem known in the art of optical flow, the flow in the x direction may be dominated by noise, while the optical flow in they direction can be clean. Knowing the uncertainty in the flow allows uncertain estimates to be down-weighted, increasing the robustness of the method.

(166) To find the variance of the optical flow, the following equation can be formulated as a posterior distribution:

(167) $P\left(v|I\right)=\exp \left(-\underset{x}{.Math.}{}_1{.Math.\frac{I}{x}{v}_x+\frac{I}{y}{v}_y+\frac{I}{t}.Math.}^2+{}_2{.Math.\frac{v}{x}.Math.}^2+{}_2{.Math.\frac{v}{y}.Math.}^2\right)$
Here, P(v|I) is a Gaussian distribution, and the mean of P(v|I) equals to the solution of the original optical flow equation for {tilde over (v)} given above. With this formulation, the variance of the optical flow (the wiggle features) can also be calculated. See the supplementary material for the detailed calculation. Let {tilde over (v)} and v be the mean and covariance, respectively, of the wiggle features computed from Equation the equation above for posterior distribution P(v|I). The L2-norm for regularization is used, in contrast to robust penalty functions such as L1 norm traditionally used by optical flow methods, since fluid objects, especially hot air or gas, do not have the clear and sharp boundaries of solid objects.

(168) Then, with the variance of the wiggle features, we can reweight the fluid flow equation can be reweighted as follows:

(169) $\tilde{u}=\arg \underset{u}{\min }\underset{x}{.Math.}{}_1{.Math.\frac{\tilde{v}}{x}{u}_x+\frac{\tilde{v}}{y}{u}_y+\frac{\tilde{v}}{t}.Math.}_{{}_v}^2+{}_2\left({.Math.\frac{u}{x}.Math.}^2+{.Math.\frac{u}{y}.Math.}^2\right)+{}_3{.Math.u.Math.}^2$
Here the squared Mahalanobis distance is denoted as x.sub..sup.2=x.sup.T.sup.1x. In this formulation, the data term is reweighted by the variance of the optical flow to robustly estimate the fluid motion: wiggle features with less certainty, such as motions measured in regions of low-contrast, or of flat or 1-dimensional structure, will have low weight in fluid flow equation. To increase the robustness, the magnitude of u can also be to avoid estimating spurious large flows. The inclusion of the above uncertainty information leads to more accurate estimation of the fluid motion, as shown in FIG. 1D(d).

(170) In practice, calculating the covariance matrix precisely for each pixel is computationally intractable, as the marginal probability distribution for each optical flow vector needs to be computed. To avoid this calculation, we concatenate all the optical flow vectors can be concatenated into a single vector and its covariance can be computed. See the supplementary material given hereinafter for the details. Also, the fluid flow equation still has a quadratic form, so the posterior distribution of the fluid flow u can be modelled as a Gaussian distribution, and its variance can be computed. This variance serves as a confidence measure in the estimated fluid motion.

(171) The goal of fluid depth estimation is to recover the depth D(x, y) of a refractive fluid object from a stereo sequence I.sub.L(x, y, t) and I.sub.R(x, y, t) (FIG. 2(b)). According to the refractive stereo constancy theorem (Theorem 2), the depth of the fluid object can be found by stereo-matching the refraction wiggles from the left and right views. So first the method discussed in the previous section can be used to calculate the mean and variance of the optical flows in the left and right views, v.sub.LN({tilde over (v)}.sub.L, .sub.L) and v.sub.RN({tilde over (v)}.sub.R, .sub.R), respectively.

(172) A discrete Markov Random Field (MRF), common for stereo matching can then be used to regularize the depth estimates. Formally, let x.sub.L, and x.sub.R be the projection of a point on the fluid object onto the left and right image plane, respectively, and define ECCV-14 submission ID 1550 9 disparity as d=x.sub.Lx.sub.R. The disparity map can be solved first by minimizing the following objective function:

(173) $\tilde{d}=\underset{d}{\min }\underset{x,y}{.Math.}f\left(v_R\left(x,y\right),v_L\left(x+d\left(x,y\right),y\right)\right)+\left({.Math.\frac{d}{x}.Math.}^2+{.Math.\frac{d}{y}.Math.}^2\right)$
where f(v.sub.R, v.sub.L) is the data term based on the observed wiggles v.sub.R and v.sub.L, and the last two terms regularize the disparity field. The inventors have found that using the L.sub.2 norm for regularization generates better results overall, better explaining the fuzzy boundaries of fluid refractive objects.

(174) As with the refractive flow, the data term can be weighted by the variance of the optical flow to make the depth estimation robust to points in a scene where the extracted wiggles are not as reliable. To achieve this, the data term, f(v.sub.R, v.sub.L), can be defined as the log of the covariance between the two optical flows from the left and right views:
f({tilde over (v)}.sub.R,{tilde over (v)}.sub.L)=log cov(v.sub.R,v.sub.L)=log|.sub.L+.sub.r|+v.sub.Rv.sub.L.sub..sub.L.sub.+.sub.R.sup.2
where v.sub.Rv.sub.L.sub..sub.L.sub.+.sub.R.sup.2=(v.sub.Rv.sub.L).sup.T(.sub.L+.sub.R).sup.1(v.sub.Rv.sub.L). This data term will assign a higher penalty to inconsistencies in the wiggle matching where the wiggles are more reliable, and a lower penalty where the wiggles are less reliable (typically where the background is less textured and the optical flow is noisy). The choice of the log of the covariance as the data term is not arbitrary, and it is the log of the conditional marginal distribution of v.sub.L and v.sub.R, given that v.sub.L and v.sub.R match. See the supplementary material given hereinafter for a detailed derivation.

(175) With calibrated cameras, the depth map, D(x, y), can be computed from the disparity map, d(x, y).

(176) Qualitative Results

(177) FIG. 7 shows several results of refractive flow analysis using a single camera.

(178) FIG. 7 contains several examples of measuring and localizing hot air radiating from several heat sources, such as people, fans, and heating vents. Refer to the supplementary materials given hereinafter for the full sequences and results.

(179) All the videos were recorded in raw format to avoid compression artifacts. To deal with small camera motions or background motions, the mean flow for each frame can be subtracted from the optical flow result. For each sequence captured, a temporal Gaussian blur was first applied to the input sequence in order to increase the SNR. The high-speed videos were captured using a high-speed camera. For some of the indoor high-speed sequences, a temporal band-stop filter was used to remove intensity variations from the light intensities modulated by AC power.

(180) The refractive flow method was first tested in controlled settings using a textured background. In the hand column, a video at 30 frames per second (fps) was taken of a person's hand right after holding a cup of hot water.

(181) Heat radiating upward from the hand was recovered. In hairdryer, a 1000 fps high-speed video of two hairdryers placed opposite to each other (the two dark shapes in the top left and bottom right are the front ends of the hairdryers) was taken, and two opposite streams of hot air flows were detected.

(182) The kettle and vents columns demonstrate the result on more natural backgrounds. In vents (700 fps), the background is very challenging for traditional background oriented schlieren (BOS) methods, as some parts of the background are very smooth or contain edges in one direction, such as the sky, the top of the dome, and the boundary of the buildings or chimneys. BOS methods rely on the motion calculated from input videos, similar to the wiggle features shown in the second row of FIG. 7, which is very noisy due to the insufficiently textured background. In contrast, the fluid flow (third row of FIG. 7) clear shows several flows of hot air coming out from heating vents. By modeling the variance of wiggles in the probabilistic refractive flow, most of the noise in the motion is suppressed. The bottom (fourth) row of FIG. 7 shows the variance of the estimated refractive flow (the square root of the determinant of the covariance matrix for each pixel).

(183) Quantitative Evaluation

(184) FIGS. 8A-8E illustrate quantitative evaluation of the refractive flow method disclosed herein. Simulated fluid density (FIG. 8A) and velocity field (FIG. 8B) were generated by Stable Fluid, a physics-based fluid flow simulation technique, and rendered on top of three different textures (FIG. 8D). The recovered velocity field from one of the simulations in which the fluid was at 320 Celsius (FIG. 8C) is similar to the ground truth velocity field (FIG. 8B). Quantitative evaluation is given in (FIG. 8E). As expected, larger temperature-related index of refraction differences between the moving fluid and the background give better flow estimates.

(185) To quantitatively evaluate the fluid velocity recovered by the disclosed refractive flow method, it on was also tested simulated sequences with precise ground truth reference. A set of realistic simulations of dynamic refracting fluid was generated using Stable Fluids, a physics-based fluid flows simulation technique, resulting in fluid densities and (2D) ground truth velocities at each pixel over time, as illustrated in FIGS. 8A-8B. The simulated densities were then used to render refraction patterns over several types of background textures with varying spatial smoothness (FIG. 8D). To render realistic refractions, it was assumed that the simulated flow was at a constant distance between the camera and background with thickness depending linearly on its density. Given a scene temperature (20 Celsius was used) and a temperature of the fluid, the amount of refraction at every point in the image can be computed. The refractive flow method previously described was then applied to the refraction sequences. The endpoint error of the recovered fluid motion at different temperature is shown in (FIG. 8E). All the sequences and results are available in the supplemental material provided hereinafter.

(186) FIGS. 9A-9D illustrate a controlled experiment to evaluate the velocities estimated using a refractive flow method. (FIG. 9A) The experiment setup. (FIG. 9B) A representative frame from the captured video. (FIG. 9C) The mean velocity of the hot air blown by the hairdryer, as computed by our method, in meters per second (m/s). (FIG. 9D) Numerical comparison of our estimated velocities with velocities measured using a velometer, for the four points marked x.sub.1-x.sub.4 in (FIG. 9C).

(187) To demonstrate further that the magnitude of the motion computed by the refractive flow method is correct, a controlled experiment as shown in FIGS. 9A-9D. A velocimeter was used to measure the velocity of hot air blown from a hairdryer, the velocities were compared with the velocities extracted by the refractive flow method. To convert the velocity on the image plane to the velocity in the real world, the hairdryer was set parallel to the imaging plane and a ruler was placed on the same plane of the hairdryer. The velocity of the flow was measured at four different locations as shown in FIGS. 9C-9D. Although the estimated flow does not match the velocimeter flow exactly, they are highly correlated (FIG. 9D) and the agreement is believed to be accurate to within the experimental uncertainties.

(188) Qualitative Results

(189) FIG. 10 illustrates fluid depth recovery from stereo sequences. The first and second rows show the input videos overlayed by wiggle features (optical flow) calculated. third line shows the disparity map (weighted by the confidence of disparity) of the fluid object.

(190) Several stereo sequences were captured using two cameras at acquiring images at 50 FPS. The two cameras were synchronized via genlock and a global shutter was used to avoid temporal misalignment. All videos were captured in 16 bit grayscale.

(191) A controlled experiment to evaluate the velocities estimated by the method. (a) The experiment setup. (b) A representative frame from the captured video. (c) The mean velocity of the hot air blown by the hairdryer, as computed using the method, in meters per second (m/s). (d) Numerical comparison of our estimated velocities with velocities measured using a velometer, for the four points marked x.sub.1-x.sub.4 in (c).

(192) The third row of FIG. 10 shows the disparity map of air flow. In the candle column are shown two plumes of hot air at different depth recovered. In lighting, three lights were positioned at different distances from the camera: the left one was the closest to the camera and the middle one was the furthest. The disparities of three plumes calculated in accordance with the method match the locations of the lights. In monitor, the method recovers the location of hot air rising from the top of a monitor. The right side of the monitor is closer to the camera, and so the hot air on the right has larger disparity. In the supplementary material given hereinafter are shown the recovered positions of the hot air as the heater rotates back and forth over time.

(193) Quantitative Evaluation

(194) A synthetic simulation was performed to demonstrate that probabilistic refractive stereo is robust to less-textured backgrounds. The simulation setup is similar to one used for the refractive flow method, and its details and the results are available in the supplementary material given hereinafter. The probabilistic framework allows for depth estimates of the flow even over image regions of only 1-dimensional structure.

(195) Supplementary Material

(196) Refractive Stereo: Quantitative Evaluation

(197) Refractive stereo results were evaluated quantitatively in two ways. First, for natural stereo sequences, recovered depths of the refractive fluid layers were compared with those of the heat sources generating them, as computed using an existing stereo method (since the depth of the actual heat sources, being solid objects, can be estimated well using existing stereo techniques). More specifically, a region on the heat source was chosen and another region of hot air right above the heat source was chosen, and the average disparities in these two regions was compared. Experiments showed that the recovered depth map of the (refractive) hot air matches well the recovered depth map of the (solid) heat source, with an average error of less than a few pixels.

(198) Second, the refractive stereo method was evaluated on simulated sequences with ground truth disparity. The simulation setup is similar to one used for the refractive flow method, except that the ground truth disparity map was manually specified as shown in FIGS. 11A-11B.

(199) To evaluate the performance of the method, four different background patterns were generated as shown in FIGS. 11A-11B: 1) weakly textured in both x and y directions, 2) strongly textured in x direction, 3) strongly textured in y direction, 4) strongly textured in both directions. All combinations of these backgrounds were generated in the left and right views to yield 44=16 stereo sequences. For comparison, a refractive stereo method that does not consider the uncertainty in the optical flow was also implemented. FIGS. 11A-11B show the disparities recovered by both refractive stereo (without uncertainty; the baseline) and probabilistic refractive stereo (with uncertainty), together with their corresponding root-mean-square error (RMSE). The simple refractive stereo worked well when the backgrounds of both the left and right views were strongly textured in both the x and y directions. When the background was weakly textured in one direction, the optical flow in that direction was noisy and thus the recovered disparity would be significantly less accurate.

(200) The probabilistic refractive stereo method was able to handle weaker textures. As long as one direction of the background was textured in both views, the disparity map was accurately recovered. For example, in the second row and second column in FIG. 1B, the optical flow in the x direction is noisy, while the optical flow in they direction is clean, because the background is very smooth in x direction. In the probabilistic refractive stereo, the optical flow results were weighted by how accurately they are estimated when calculating the disparity. Therefore, the optical flow in the x direction would have a smaller effect on the final result than the optical flow in they direction, and the disparity map would be calculated correctly. This demonstrates the robustness of the method to partially textured backgrounds.

(201) FIG. 12 is similar to FIG. 6B but also illustrates quantities relevant to a proof of refractive constancy. As described above, the refractive wiggle is defined as:

(202) Definition 1 (Refraction wiggle) Let x be a point on the refractive fluid layer, x be the intersection between the background and the light ray passing through x and the center of projection at time t, and t be a short time interval. Then the wiggle of x from time t to t+t is the shift of the projection of x on the image plane during this time.

(203) This definition leads to the following two refractive constancy theorems.

(204) Theorem 1 (Refractive flow constancy) Suppose the fluid object does not change its shape and index of refraction during a short time interval [t.sub.1, t.sub.2]. Then for any point on the fluid object, its wiggle v(t.sub.1) at t.sub.1 equals its wiggle v(t.sub.2) at t.sub.2.

(205) Proof Let x.sub.t.sub.1 be the location of a point on the fluid object at time t.sub.1, and x.sub.t.sub.2 be the location of the same point at time t.sub.2. It will be shown that the wiggles at these two points and times are the same.

(206) Because wiggles are defined by shifts in the image plane, rays are first traced to determine what points in the image plane correspond to these locations on the fluid object. At time ti, i=1, 2, an undistorted ray is emitted from the center of the projection to point x.sub.t.sub.i with the angle .sub.t.sub.i with respect to vertical (solid lines in FIG. 2). This ray is bent by refraction as it passes through the fluid and travels at an angle .sub.t.sub.i with respect to the vertical before finally intersecting the background at x.sub.t.sub.i.

(207) At a successive time t.sub.i+t, the fluid object moves to a new location (dashed gray blob in FIG. 12). The background point x.sub.t.sub.i then corresponds to a different point x.sub.t.sub.i+t on the fluid object and a different point x.sub.t.sub.i+t on the image plane. The ray from the background to the fluid object now makes an angle of .sub.t.sub.i+t with the vertical and the same ray as it travels from the fluid object to the center of projection makes an angle of .sub.t.sub.i+t with the vertical. The wiggle at time t.sub.i is the distance between x.sub.t.sub.i and x.sub.t.sub.i+t, denoted as {right arrow over (x.sub.t.sub.1x.sub.t.sub.1+t)} (blue and red arrows in FIG. 12). It can be proven that the wiggles at time t.sub.1 and t.sub.2 are equal to each other. That is,
{right arrow over (x.sub.t.sub.1x.sub.t.sub.1+t)}={right arrow over (x.sub.t.sub.2x.sub.t.sub.2+t)}

(208) To prove this, it will first be shown that

(209) $\stackrel{x_{t_1+t}^{}}{.fwdarw.}\mathrm{and}\stackrel{x_{t_2+t}^{}}{.fwdarw.}$
same and are the same point on the fluid object, or equivalently, the shifts x.sub.t.sub.ix.sub.t.sub.1+t and x.sub.t.sub.2x.sub.t.sub.2+t are equal. By the definition of wiggle, it is known that x.sub.t.sub.1 and x.sub.t.sub.1 refer to the same point on the refractive object. It can be proven that x.sub.t.sub.1+t and x.sub.t.sub.2+t are equal.

(210) Let z, z, and z, respectively, be the depths of the image plane, the fluid layer, and the background from the center of projection (FIG. 2). To find the relationship between x.sub.t.sub.1+t and x.sub.t.sub.2+t, x.sub.t.sub.2+t is solved for in terms of x.sub.t.sub.1, x.sub.t.sub.2, and x.sub.t.sub.1+t. A simple geometric calculation shows (assuming that all the incident angle t.sub.i is small, so that a.sub.t.sub.isin(.sub.t.sub.i)):

(211) $\frac{x_{t_i}^{}-o}{z^{}}={}_{t_i}^{}\left(i=1,2\right),\frac{x_{t_i}^{}-x_{t_i}^{}}{z^{}-z^{}}={}_{t_i}^{}\left(i=1,2\right).$

(212) The relationship between .sub.t.sub.i and .sub.t.sub.i is due to refraction and is determined by Snell's law. The first order paraxial approximation for small incident angles, as commonly used in the optics literature. The result is that Snell's refraction law is approximately an additive constant independent of the original ray direction.
.sub.t.sub.i.sub.t.sub.i+.sub.t(x.sub.t.sub.i).

(213) The angle difference t(x.sub.t.sub.i) depends only on the normal direction of fluid object's surface at x.sub.t.sub.i and the fluid object's index of refraction, and is therefore independent of the incoming ray direction .sub.t.sub.i. Because it is assumed that the fluid object does not change shape and index of refraction in the short time interval, the difference in angles also remains constant during the short time interval. Subtracting the first equation in this section from the second one results in:

(214) ${}_t={}_{t_i}^{}-{}_{t_i}^{}=\frac{x_{t_i}^{}-x_{t_i}^{}}{z^{}-z^{}}-\frac{x_{t_i}^{}-o}{z^{}}\left(i=1,2\right),$

(215) From which for x.sub.t.sub.i can be solved for in terms of the other variables:

(216) $x_{t_i}^{}-o=\frac{z^{}}{z^{}}\left(x_{t_i}^{}-o\right)-\frac{z^{}\left(z^{}-z^{}\right)}{z^{}}{}_t\left(x_{t_i}^{}\right)\left(i=1,2\right).$

(217) Similarly, we can solve for x.sub.t.sub.i+t:

(218) 0$x_{t_i}^{}-o=\frac{z^{}}{z^{}}\left(x_{t_i}^{}-o\right)-\frac{z^{}\left(z^{}-z^{}\right)}{z^{}}{}_t\left(x_{t_i}^{}\right)\left(i=1,2\right).$

(219) Subtracting the previous two equations from each other yields a difference equation:

(220) $x_{t_i+t|}^{}-x_t^{}=-\frac{z^{}\left(z^{}-z^{}\right)}{z^{}}\left({}_{t_i+t}\left(x_{t_i+t}^{}\right)-{}_{t_i}\left(x_{t_i}^{}\right)\right)\left(i=1,2\right).$

(221) Since x.sub.t.sub.1 and x.sub.t.sub.2 are the same point on the fluid object, differences between incidence angles and refraction angles should be same at x.sub.t.sub.1 and x.sub.t.sub.2, that is .sub.t.sub.i(x.sub.t.sub.1)=.sub.t.sub.2(x.sub.t.sub.2). Therefore, subtracting the previous equation with i=2 from the previous equation with i=1, yields

(222) $\left(x_{t_1+t}^{}-x_{t_1}^{}\right)-\left(x_{t_2+t}^{}-x_{t_2}^{}\right)=-\frac{z^{}\left(z^{}-z^{}\right)}{z^{}}\left({}_{t_1+t}\left(x_{t_1+t}^{}\right)-{}_{t_2+t}\left(x_{t_2+t}^{}\right)\right).$

(223) The quantity x.sub.t.sub.2+t from the previous equation can then be solved for. If x.sub.t.sub.2+t and x.sub.t.sub.2+t are also the same point on the fluid object, then the left hand side (LHS) of the previous equation is zero as x.sub.t.sub.1+tx.sub.t.sub.1=x.sub.t.sub.2+tx.sub.t.sub.2, and the right hand side (RHS) of the previous equation is also zero as .sub.t.sub.1+t(x.sub.t.sub.1+t)=.sub.t.sub.2+t(x.sub.t.sub.2+t). This shows that x.sub.t.sub.2+t=x.sub.t.sub.1+tx.sub.t.sub.1+x.sub.t.sub.2 is the solution to the difference equation above, and therefore x.sub.t.sub.1+t and x.sub.t.sub.2+t are also the same point on the fluid object.

(224) Finally, it is proven that wiggles at t.sub.1 and t.sub.2 are equal. By the similar triangle formula, we have:

(225) $\stackrel{.Math.}{x_{t_1}{x}_{t_1+t}}=\frac{z}{z^{}}\stackrel{.Math.}{x_{t_1}^{}{x}_{t_1+t}^{}}=\frac{z}{z^{}}\stackrel{.Math.}{x_{t_2}^{}{x}_{t_2+t}^{}}=\stackrel{.Math.}{x_{t_2}{x}_{t_2+t}}.$

(226) Theorem 2 (Refractive stereo constancy) Suppose there are n2 cameras imaging a refractive fluid object, and they are all parallel and close to each other. Then at any time t, and for any point on the fluid object, the corresponding wiggle in all the cameras is the same.

(227) FIG. 13 is similar to FIG. 6C but also illustrates quantities relevant to a proof of stereo refractive constancy.

(228) Proof Similar to the previous proof, let xt be a point on the fluid object. Tracing rays, at time t, an undistorted ray is emitted from the center of the projection o.sub.j(j=1, 2) to the point point x.sub.t. It is refracted by the fluid object and intersects the background at x.sub.j. At a successive time t+t, the fluid object moves, and the light ray from the points on the background to the center of the projection now goes through points x.sub.t+t,j on the fluid object and x.sub.t.sub.i+t,j on the image plane. It will be shown that the wiggle features

(229) $\stackrel{.Math.}{x_{t,1}{x}_{t+t,1}}\mathrm{and}\stackrel{.Math.}{x_{t,2}{x}_{t+t,2}}$
are equal in both views.

(230) As in the previous proof, it is first shown, that x.sub.t+t,1=x.sub.t+t,2. Following a similar derivation as in the previous proof results in:

(231) $x_{t+t,j}^{}-x_t^{}=-\frac{z^{}\left(z^{}-z^{}\right)}{z^{}}\left({}_{t+t}\left(x_{t+t,j}^{}\right)-{}_t\left(x_t^{}\right)\right).$

(232) Then, subtracting the previous equation with j=2 from the previous equation with j=1, we have:

(233) $x_{t+t,1}^{}-x_{t+t,2}^{}=-\frac{z^{}\left(z^{}-z^{}\right)}{z^{}}\left({}_{t+t}\left(x_{t+t,1}^{}\right)-{}_{t+t}\left(x_{t+t,2}^{}\right)\right).$

(234) By the same logic as in the previous proof, when x.sub.t+t,1=x.sub.t+t,2, both the LHS and RHS of the previous equation are equal to 0. Therefore, x.sub.t+t,1=x.sub.t+t,2 is the solution to the previous equation.

(235) Thus, x.sub.t+t,1 and x.sub.t+t,2 are the same point.

(236) Finally, it can be proven that wiggles from two views are equal:

(237) $\stackrel{.Math.}{x_{t+t,1}{x}_{t,1}}=\frac{z}{z^{}}\stackrel{.Math.}{x_{t+t,1}^{}{x}_t^{}}=\frac{z}{z^{}}\stackrel{.Math.}{x_{t+t,2}^{}{x}_t^{}}=\stackrel{.Math.}{x_{t+t,2}{x}_{t+t,2}}.$
Calculating Fluid Flow Efficiently

(238) It was previously shown that the probabilistic refractive flow method consists of two steps. First, the mean {tilde over (v)} and the variance .sub.v of the wiggle features v are solved for from the following Gaussian distribution:

(239) $\stackrel{.Math.}{x_{t+t,1}{x}_{t,1}}=\frac{z}{z^{}}\stackrel{.Math.}{x_{t+t,1}^{}{x}_t^{}}=\frac{z}{z^{}}\stackrel{.Math.}{x_{t+t,2}^{}{x}_t^{}}=\stackrel{.Math.}{x_{t+t,2}{x}_{t+t,2}}.$

(240) To solve for the mean and variance of flow from the previous equation, let the V be the vector formed by concatenating all the optical flow vectors in one frame. That is, V=( . . . , v(x), . . . ). Also, let us represent the previous equation in information form P(v|I)=exp(V.sup.TJV+h.sup.TV), where h and J can be calculated from the previous equation. Then the mean of V is {tilde over (v)}=J.sup.1h and covariance of V is =J.sup..

(241) In the second step, the fluid flow is calculated by minimizing the following optimization problem based on the mean and variance of the wiggle features computed in the first step.

(242) $\tilde{u}=\arg \underset{u}{\min }\underset{x}{\overset{}{.Math.}}{}_1{.Math.\frac{\tilde{v}}{x}{u}_x+\frac{\tilde{v}}{y}{u}_y+\frac{\tilde{v}}{t}.Math.}_{{}_v}^2+{}_2\left({.Math.\frac{u}{x}.Math.}^2+{.Math.\frac{u}{y}.Math.}^2\right)+{}_3{.Math.u.Math.}^2$

(243) Calculating the covariance of each wiggle feature requires inverting the information matrix J. This step will be slow if the matrix is large. To avoid this time consuming inversion, we make a slight change to the fluid flow objective function. Let {tilde over (V)}.sub.x, {tilde over (V)}.sub.y, and {tilde over (V)}.sub.t be the vectors formed by concatenating all the partial derivatives of mean wiggle features in a frame, that

(244) 0${\tilde{V}}_x=(\mathrm{.Math.},\frac{v}{x}\left(x\right),\mathrm{.Math.}),{\tilde{V}}_y=(\mathrm{.Math.},\frac{v}{y}\left(x\right),\mathrm{.Math.}),\mathrm{and}$${\tilde{V}}_t=(\mathrm{.Math.},\frac{v}{t}\left(x\right),\mathrm{.Math.}).$
Similarly, let U.sub.x, U.sub.y be the vectors formed by concatenating all the x-components and y-components of u in a frame respectively. Then the refractive flow can be calculated as follows:

(245) $\underset{U}{\min }{{}_1\left({\tilde{V}}_x.Math.U_x+{\tilde{V}}_y.Math.U_y+{\tilde{V}}_t\right)}^{T}J\left({\tilde{V}}_z.Math.U_x+{\tilde{V}}_y.Math.U_y+{\tilde{V}}_t\right)+{}_2\left({.Math.D_x{U}_x.Math.}^2+{.Math.D_y{U}_x.Math.}^2+{.Math.D_x{U}_y.Math.}^2+{.Math.D_y{U}_y.Math.}^2\right)+{}_3{.Math.U.Math.}^2$

(246) where D.sub.x and D.sub.y are the partial derivative matrices to x and y respectively. The smoothness term of the previous equation is exactly the same as that in the equation for above, and the data term of the equation for above is
(V.sub.x.Math.U.sub.x+V.sub.y.Math.U.sub.y+V.sub.t).sup.TJ(V.sub.x.Math.U.sub.x+V.sub.y.Math.U.sub.y+V.sub.t)=V.sub.x.Math.U.sub.x+V.sub.y.Math.U.sub.y+V.sub.t.sub.J-1=V.sub.x.Math.U.sub.x+V.sub.y.Math.U.sub.y+V.sub.t.sub.,
which is also similar to the data term in the equation for above except that it jointly considers all the wiggle vectors in a frame. Therefore, this change will not affect the result too much, but the method is more computationally efficient, as we never need to compute J.sup.1. The term never appears in the refractive flow expression above.
Probabilistic Interpretation in the Refractive Stereo

(247) In this section, it is shown that the data term defined in Section 5 for refractive stereo is equal to the negative log of the conditional marginal distribution. Let v.sub.R(x)N(v.sub.R(x),.sub.R(x)) and v.sub.L(x+d(x))N(v.sub.L(x+d(x)), .sub.L(x+d(x))) be the optical flow from the left and right views, where v.sub.L and v.sub.R are the means of the optical flow and and .sub.L and .sub.R are the variances of optical flow. Then the data term of the refractive stereo is

(248) $f\left(v_R,v_L\right)=-\log \mathrm{cov}\left(v_L,v_R\right)=-\log {}_v{P}_L\left(v\right)P_R\left(v\right)\mathrm{dv}=-\log {}_{v}N\left(v;{\stackrel{\_}{v}}_L,{}_L\right)N\left(v;{\stackrel{\_}{v}}_R,{}_R\right)\mathrm{dv}=\frac{1}{2}\log .Math.{}_L+{}_R.Math.+\frac{1}{2}{.Math.{\stackrel{\_}{v}}_R-{\stackrel{\_}{v}}_L.Math.}_{{}_L+{}_R}^2+\mathrm{const}$

(249) Where N(v; v.sub.L, .sub.L) is shorthand notation for the Gaussian distribution probability density function

(250) $N\left(v;\stackrel{\_}{v},\right)=\frac{1}{2\sqrt{.Math..Math.}}\exp \left(-\frac{1}{2}{\stackrel{\_}{v}}^T{}^{-1}\stackrel{\_}{v}\right)$

(251) Recall that the optical flow calculated by the method is degraded by noise. Specifically, let v(x) be the ground truth optical flow from the right view at x. The mean optical flow from the right view (or left view) calculated by the method equals the ground truth optical flow plus Gaussian noise with zero-mean and variance equal to .sub.R(or .sub.L), that is:
P(v.sub.R(x)|v(x))=N(v.sub.R(x);v(x),.sub.R(x)),P(v.sub.L(x+d(x))|v(x),d(x))=N(v.sub.L(x+d(x));v(x),.sub.L(x+d(x)))

(252) To evaluate the probability of d, let us consider the marginal distribution

(253) $P\left({\stackrel{\_}{v}}_R\left(x\right),{\stackrel{\_}{v}}_L\left(x+d\left(x\right)\right)|d\left(x\right)\right)={}_{v}P\left({\stackrel{\_}{v}}_R\left(x\right),{\stackrel{\_}{v}}_L\left(x+d\left(x\right)\right),v\left(x\right)|d\left(x\right)\right)\mathrm{dv}\left(x\right)={}_{v}P\left({\stackrel{\_}{v}}_R\left(x\right)|v\left(x\right)\right)P\left({\stackrel{\_}{v}}_L\left(x+d\left(x\right)\right)|v\left(x\right),d\left(x\right)\right)P\left(v\right)\mathrm{dv}\left(x\right)={}_{v}N\left({\stackrel{\_}{v}}_R\left(x\right);v\left(x\right),{}_R\left(x\right)\right)N\left({\stackrel{\_}{v}}_L\left(x+d\left(x\right)\right);v\left(x\right),{}_L\left(x+d\left(x\right)\right)\right)P\left(v\right)\mathrm{dv}\left(x\right)$

(254) Assuming that P(v) has an uniform prior, we have:

(255) $\log P\left({\stackrel{\_}{v}}_R\left(x\right),{\stackrel{\_}{v}}_L\left(x+d\left(x\right)\right)|d\left(x\right)\right)=\log {}_{v}N\left({\stackrel{\_}{v}}_R\left(x\right);v,{}_R\left(x\right)\right)N\left({\stackrel{\_}{v}}_L\left(x+d\left(x\right)\right);v,{}_L\left(x+d\left(x\right)\right)\right)P\left(v\right)dv=\log {}_{v}N\left({\stackrel{\_}{v}}_R\left(x\right);v,{}_R\left(x\right)\right)N\left({\stackrel{\_}{v}}_L\left(x+d\left(x\right)\right);v,{}_L\left(x+d\left(x\right)\right)\right)dv+\mathrm{const}=\log {}_{v}N\left(v;{\stackrel{\_}{v}}_R\left(x\right),{}_R\left(x\right)\right)N\left(v;{\stackrel{\_}{v}}_L\left(x+d\left(x\right)\right),{}_L\left(x+d\left(x\right)\right)\right)dv+\mathrm{const}=-f\left(v_R,v_L\right)+\mathrm{const}$

(256) Therefore, the data term is equal to the negative log of conditional marginal distribution (plus a constant).

(257) While this invention has been particularly shown and described with references to example embodiments thereof, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the scope of the invention encompassed by the appended claims.