SYSTEM AND METHOD FOR REPRESENTING THE TACTILE INTERACTION EMPLOYED BY AN ARRAY OF ULTRASOUND TRANSDUCERS

Abstract

A system and method for rendering tactile interaction using an ultrasound transducer array by receiving a target field; computing one or more paths of pressure focal points of varying intensity such that the difference between a result field produced by the one or more paths of pressure focal points and the target field is minimized; and controlling an ultrasound transducer array in spatiotemporal modulation mode to simultaneously render the one or more computed paths of pressure focal points. The system and method for rendering tactile interaction using an ultrasound transducer array accurately represents dynamic interactions produced in a virtual environment with a varying pressure or force field.

Claims

1. A method for rendering tactile interaction using an ultrasound transducer array comprising: receiving a target field, wherein the target field is a force field or a pressure field; computing one or more paths of pressure focal points of varying intensity such that the difference between a result field produced by the one or more paths of pressure focal points and the target field is minimized; and controlling an ultrasound transducer array in spatiotemporal modulation mode to simultaneously render the one or more computed paths of pressure focal points.

2. The method of claim 1, wherein the step of computing one or more paths of pressure focal points comprises: obtaining one or more initial paths of pressure focal points; and optimizing the one or more initial paths of pressure focal points by minimizing a set of cost terms based on the target field and the result field produced by the one or more paths of pressure focal points.

3. The method of claim 2, wherein the set of cost terms is minimized by iterating steps of gradient descent, wherein the set of cost terms includes an intensity cost term, c.sub.p, to maximize field intensity.

4. (canceled)

5. The method of claim 2, wherein the set of cost terms includes one of the following: a length cost term, c.sub.1 to comply with a target path length, L, an intersection cost term, c.sub.i, to prevent path samples from getting closer to a fall-off distance σ, or a bending cost term, c.sub.b, to obtain low-curvature paths.

6. (canceled)

7. (canceled)

8. The method of claim 2, wherein the step of obtaining one or more initial paths of pressure focal points comprises: obtaining, from the target field, a plurality of target points with an associated position and target value; and determining the one or more initial paths of pressure focal points by selecting one or more sequences of target points.

9. The method of claim 1, wherein the step of computing one or more paths of pressure focal points comprises: obtaining, from the target field, a set of sample points with an associated target value; and determining one or more initial paths of pressure focal points by: computing the shortest closed path that visits all the points of the set of sample points; and recursively splitting the set of sample points into subsets and computing the shortest closed path for each subset until the path length satisfy a maximum path length constraint for all the subsets.

10. The method of claim 9, wherein the set of sample points are obtained by applying a cluster analysis on a plurality of target points of the target field.

11. The method of claim 10, wherein the cluster analysis is a weighted k-means algorithm, wherein the result field is obtained by applying a spatial smoothing function to the intensity of the focal points, and wherein the spatial smoothing function is defined as: $p\left(x\right)=\frac{1}{\gamma }{p}_i\varphi \left(.Math.x-x_i.Math.\right)$ wherein p(x) is the intensity of the result field at position x, x.sub.i is the closest position to x in a path, p.sub.i is the intensity of the target field at position x.sub.i, ϕ is a Gaussian function and γ is a heuristic gain.

12. (canceled)

13. (canceled)

14. The method of claim 1, further comprising projecting the target field onto a plane to obtain a plurality of target points defined by a 2D position (x.sub.i) and a target value, undoing the projection to obtain the one or more paths of pressure focal points, refining the target field by computing the intersection of the received target field with a surface of a body or a part of a body of a user, tracking the position of the body or the part of the body of the user.

15. (canceled)

16. (canceled)

17. (canceled)

18. The method of claim 1, wherein the target field is generated from an interaction of a user with a virtual environment.

19. The method of claim 1, wherein the target field is generated from an interaction of a user with a real environment remote from the user location.

20. The method of claim 19, wherein each path of pressure focal points is a 3D curve defining a closed path.

21. (canceled)

22. The method of claim 1, wherein each path of pressure focal points satisfies a maximum path length constraint.

23. (canceled)

24. The method of claim 1, further comprising computing an interactive simulation and extracting the target field from the simulation.

25. A system for rendering tactile interaction using an ultrasound transducer array, wherein the system comprises a data processing device configured to: receive a target field, wherein the target field is a force field or a pressure field; compute one or more paths of pressure focal points of varying intensity such that the difference between a result field produced by the one or more paths of pressure focal points and the target field is minimized; and control an ultrasound transducer array in spatiotemporal modulation mode to simultaneously render the one or more computed paths of pressure focal points.

26. The system of claim 25, further comprising a driver circuit for actuating the ultrasound transducer array.

27. (canceled)

28. The system of claim 25, further comprising a tracking device for tracking the position of the body or the part of the body of the user.

29. (canceled)

30. The system of claim 25, wherein the target field is generated from an interaction of a user with a virtual environment.

31. The system of claim 25, wherein the data processing device is configured to compute an interactive simulation and extract the target field from the simulation.

32. (canceled)

33. A computer program product for rendering tactile interaction using an ultrasound transducer array comprising computer usable program code for, when executed on a processor, performing the steps of the method defined in claim 1.

34. (canceled)

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0034] A series of drawings which aid in better understanding the invention and which are expressly related with an embodiment of said invention, presented as a non-limiting example thereof, are very briefly described below.

[0035] FIG. 1 depicts a flow diagram of a method of rendering tactile interaction using an ultrasound transducer array according to an embodiment of the present invention.

[0036] FIG. 2 shows the steps of computing paths of pressure focal points according to an embodiment.

[0037] FIG. 3 shows the steps of computing paths of pressure focal points according to another embodiment.

[0038] FIG. 4 is a flow diagram for the step of obtaining initial paths of pressure focal points in FIG. 2, according to an embodiment.

[0039] FIG. 5A shows the target pressure points of a target field produced on the user's hand and extracted from a fluid simulation. FIG. 5B represents the down-sampled pressure points (clusters) of the target field. FIG. 5C shows the optimal initial path for the full set of pressure sample points. FIG. 5D shows the optimal paths after splitting the set of sample points into subsets to satisfy a maximum path length constraint. FIG. 5E depicts the refinement of the paths to maximize pressure intensity. FIG. 5F shows the resulting reconstructed pressure field obtained by the optimization algorithm.

[0040] FIG. 6 depicts a system for rendering tactile interaction using an ultrasound transducer array according to an exemplary embodiment.

[0041] FIG. 7 depicts the virtual representation of a user's hand interacting with virtual plumes of smoke used in a dynamic fluid simulation within a virtual environment.

[0042] FIG. 8A depicts the pressure target field produced by the smoke plumes on the surface of the virtual hand. FIG. 8B illustrates the reconstructed pressure field produced by STM rendering of dynamically optimized focal point paths according to the present method.

[0043] FIGS. 9 and 10 depict experimental virtual scenarios (with one or four plumes of smoke, respectively) and representative examples of target and reconstructed fields using the STM-based method of the present invention and an AM-based method of the prior art.

DESCRIPTION OF A PREFERRED EMBODIMENT OF THE INVENTION

[0044] When using a spatiotemporal modulation (STM) in the activation of the transducers of an ultrasound transducer array, a focal point traverses a path in space. The path P can be formally described as a time-dependent position: P=x(t)∈ .sup.3. The perceived intensity of a focal point path is maximized under constructive interference between the motion of the focal point and the propagation of skin waves (Frier et al. [6]). This happens for focal point speeds between 5 m/s and 8 m/s; therefore, a reference speed v within said margin is selected for the rendering algorithm. In an embodiment, the selected reference speed v is 7 m/s. Furthermore, to ensure constructive interference on the complete path, the paths are preferably designed as closed paths (i.e., the path P is a closed 3D curve).

[0045] Frier et al. [6] assumes that the length of the path is given; therefore, the frequency at which the path is repeated cannot be independently controlled, and depends on the traversal speed and the length of the path. According to preliminary experiments, it has been observed that this is acceptable up to a maximum path length. Beyond that length, the frequency at which the path is repeated is too low, and the stimulus is no longer perceived as continuous. To determine the minimum acceptable frequency (i.e., the maximum acceptable path length), circles of different radii at the reference traversal speed of 7 m/s were rendered, finding that a minimum frequency of 50 Hz (i.e., a maximum path length L=140 mm) is a safe bound to ensure that the stimulus is perceived as continuous.

[0046] Normally, the haptic devices (i.e., the ultrasound transducer arrays) allow STM of multiple focal points simultaneously. Each focal point can traverse a different path, with all focal points traveling at the reference speed v, and all paths satisfying the maximum length constraint L. From preliminary experiments, the number of simultaneous focal points can be limited to a certain number of simultaneous paths (in an embodiment, to four paths). More focal points may reach larger coverage, but at the price of notable degradation of perceived intensity.

[0047] In previous works STM is used to render 3D curves, wherein the intensity of the radiation pressure of the moving focal point is kept constant along such curves. In the present invention a temporally and spatially varying pressure field is rendered, wherein the intensity of the radiation pressure along the path adapts locally to the intensity of the pressure field. A path is characterized with a position-dependent pressure intensity p(x).

[0048] As a focal point cycles multiple times through the same position x.sub.i, the rendered pressure intensity p(x.sub.i)=p.sub.i is the same on all cycles. If the path is repeated frequently enough (e.g., at more than 50 Hz), the rendered radiation pressure produces a persistent tactile perception. This is equivalent to applying a time-invariant pressure at each position along the path, with its effective magnitude a fraction of the rendered pressure. With this assumption, during a time window the effective pressure can be considered to be a spatially varying but temporally invariant field, i.e., a quasi-static pressure field.

[0049] The effective pressure of STM rendering is not a simple time-average of the rendered pressure. According to an embodiment, a perceptual heuristic is followed to approximate this magnitude: rendering the same stimuli using amplitude modulation (AM) and STM, and asking subjects to tune the gain γ of STM until the peak perceptual intensity is similar. In practice, a gain γ=1.4 is used.

[0050] Focal points exhibit a smooth fall-off determined by the wavelength of the ultrasound signal (e.g., 8.6 mm for the 40 kHz of a test ultrasound device). As shown by Hoshi et al. [2], this fall-off can be approximated well by a Gaussian function ϕ. Based on this finding, together with the heuristic gain γ, the effective quasi-static pressure field p(x) produced by a focal point path can be approximated as:

[00002]$\begin{array}{cc}p\left(x\right)=\frac{1}{\gamma }{p}_i\varphi \left(.Math.x-x_i.Math.\right)=\frac{1}{\gamma }{p}_i{e}^{-\frac{{.Math.x-x_i.Math.}^2}{2{\sigma }^2}};& \left(1\right)\end{array}$

where x.sub.i is the closest position to x in the path. The standard deviation σ of the Gaussian fall-off is set to the same value as the wavelength of the ultrasound signal (e.g., 8.6 mm for 40 Khz).

[0051] The quasi-static pressure field assumption simplifies the rendering algorithm of the present invention. Given a target pressure field, obtained for instance from a fluid simulation, the focal point paths whose quasi-static pressure field best reconstructs the target field are searched for. To ensure that the quasi-static pressure field assumption is valid, this search preferably fulfills two constraints, namely that each focal point travels at the reference speed v and the length of each path is equal or shorter than the maximum path length L.

[0052] Given a target pressure field or a target force field, focal point paths that produce a best-matching quasi-static pressure field are sought by a path routing optimization algorithm. FIG. 1 depicts, according to an embodiment of the present invention, a flow diagram of a method 100 of rendering tactile interaction using an ultrasound transducer array that implements a path routing optimization algorithm.

[0053] The method 100 comprises receiving 110 a target field 102. The target field, which is either a force field or a pressure field, may be computed externally to the optimization method (e.g. by an external simulation). The method also comprises computing 120 one or more paths of pressure focal points of varying intensity such that the difference between a result field produced by the one or more paths of pressure focal points and the target field is minimized. Finally, the method comprises controlling 130 an ultrasound transducer array in spatiotemporal modulation mode to simultaneously render the one or more computed paths of pressure focal points.

[0054] In an embodiment, depicted in FIG. 2, the step of computing 120 one or more paths of pressure focal points comprises the following steps: obtaining 202 one or more initial paths of pressure focal points; and optimizing 204 the one or more initial paths of pressure focal points by minimizing a set of cost terms based on the target field and the result field produced by the one or more paths of pressure focal points. Alternatively or in addition to the steps of the embodiment of FIG. 2, the step of computing 120 one or more paths of pressure focal points may comprise the steps depicted in the flow diagram of FIG. 3, obtaining 302, from the target field, a set of sample points with an associated target value and determining 308 one or more initial paths of pressure focal points using these sample points. In an embodiment, the set of sample points are obtained by applying a cluster analysis on the target field 102 (in particular, a cluster analysis on a plurality of target points of the target field). The cluster analysis may be a weighted k-means algorithm, wherein the weights considered are the target value of the target points. The one or more initial paths of pressure focal points are determined by computing 304 the shortest closed path that visits all the points of the set of sample points to obtain an initial path and, if the length of the initial path exceeds a maximum path length L, recursively splitting 306 the set of sample points into subsets and computing the shortest closed path for each subset until the path length of all subsets satisfy a maximum path length constraint.

[0055] If the embodiments of FIGS. 2 and 3 are used in combination, the step of obtaining 202 one or more initial paths of pressure focal points of FIG. 2 is implemented according to the steps shown in FIG. 3. Alternatively, as illustrated in the flow diagram of FIG. 4, the step in FIG. 2 of obtaining 202 one or more initial paths of pressure focal points may comprise obtaining 402, from the target field, a plurality of target points with an associated position and target value; and determining 404 the one or more initial paths of pressure focal points by selecting one or more sequences of target points. This way, the paths may be initialized according to different initialization methods, such as randomly (e.g., FIG. 4), or applying clustering and shortest path calculation (FIG. 3). A good initialization makes the optimization algorithm more efficient and/or converging to a better solution. The paths are calculated by optimization, and there are optimization methods that need an initialization and then apply a refinement, and optimization methods performing sampling of paths and selecting the best path.

[0056] According to an embodiment, the optimization problem can be solved in two steps, at two different resolutions, searching for paths that maximize coverage and integrated pressure intensity subject to the path length constraint. The first step is a coarse step performing a global search (path initialization), while the second step is a fine step performing a local search (path refinement). In another embodiment, only one of these two steps may be necessary. For example, an embodiment may include using just the first step if the accuracy provided is sufficient, or using only the second step if the computational power allows running the fine step directly.

[0057] An interactive simulation may be computed to extract a target pressure field (or a target force field). For instance, an interactive fluid simulation can be used [8], tracking the user's hand and modelling the hand as a moving obstacle in a 3D simulation of a gaseous medium. As described in [8], fluid dynamics can be modelled using incompressible Euler equations discretized on a 3D Eulerian grid, with semi-Lagrangian advection and massively parallel Jacobi relaxation for the pressure solve. The fluid simulation can be executed on a GPU for maximum performance.

[0058] The method may comprise refining the target field by computing the intersection of the received target field with a surface of a body or a part of a body of a user. To that end, the method may comprise tracking the position of the body or the part of the body of the user. For example, to define the target pressure field the hand may be voxelized and voxel positions are selected such that they are visible from the side of the domain that corresponds to the location of the ultrasound device. To simplify the path optimization problem a plane may be fit to the voxel positions, so that the voxel positions are projected onto the plane, making path optimization a 2D problem. Formally, the target pressure field can be described by a set of pressure target points T defined by 2D positions x.sub.i and their corresponding target pressure values p*(x.sub.i):

T={(x.sub.i ∈ .sup.2, p*(x.sub.i))}.

[0059] The first coarse step includes the path initialization. To ensure high computational performance, the initialization of the paths preferably uses only a representative subset of the pressure target points (although the whole set of pressure target points T may also be considered, requiring higher computational resources). FIGS. 5A to 5D depicts the initialization step according to the embodiment of FIG. 3. Starting with a set (or subset) of pressure target points 502 of the target field 102 (e.g., pressure target points on the user's hand extracted from a fluid simulation, as shown in FIG. 5A), a set 504 of pressure sample points 506 with an associated target value are obtained 302 (FIG. 5B) by applying weighted k-means clustering to the set of pressure target points 502. The clustering analysis produces a down-sampled target pressure field formed by the pressure sample points 506.

[0060] Given the set 504 of pressure sample points 506, a set of closed paths that visit all the points is sought, subject to the maximum path length L. The optimal solution to this problem may require an arbitrarily large number of paths; however, as previously noted, it is advisable to limit the number of paths to a certain number of paths (e.g., four paths). Consequently, the resulting paths may fail to visit all the pressure sample points, and an optimal set must be selected.

[0061] This problem is solved iteratively. An optimal path (shortest closed path 508, FIG. 5C) that visits all the pressure sample points 504 is first computed 304. If the shortest closed path 508 is too long, the set 504 of pressure sample points 506 are split into two subsets and optimal paths (shortest closed paths 508 visiting all the points of the subset 510) are separately computed. The subsets 510 of pressure sample points 506 are recursively split 306 until all shortest closed paths of all the subsets 510 satisfy the maximum length constraint L. If the number of resulting paths is larger than four, the four paths with highest integrated pressure are kept (as in the example of FIG. 5D), thereby obtaining the initial paths 512.

[0062] The operations to compute an optimal path for a set of points and to split a set of points into further subsets are hereinafter described according to an embodiment, wherein the set of points are either a set 504 or a subset 510 of pressure sample points 506.

[0063] Given a set of points, finding the shortest closed path 508 that visits all the points corresponds to the traveling salesman problem. This problem can be solved using the 2-opt algorithm [9], which admits closed paths. The computational cost of 2-opt sets an upper bound on the size of the set 504 of pressure sample points 506. In an embodiment, a maximum of 50 points is set. Therefore, the weighted k-means clustering step is run with 50 or fewer clusters.

[0064] To split a set of points, the direction of maximum spread is found, the points along this direction are bounded, and a splitting plane orthogonal to the direction is placed at the midpoint of the two bounds. To find the direction of maximum spread, the covariance matrix of the points, weighted by their pressure value, is computed. The direction of maximum spread corresponds to the eigenvector with highest eigenvalue.

[0065] The second fine step includes the path refinement. After initialization, the paths pass through pressure clusters and fulfill the maximum length constraint. However, due to their coarse sampling, the initial paths are not optimally aligned with pressure peaks and ridges. Path refinement is executed at a higher resolution, by up-sampling each path to N points. In an embodiment, N is set to 20, which sets points 7 mm apart from each other, i.e., the distance traveled by a focal point in 1 ms.

[0066] During refinement, the goal is to maximize pressure intensity by upsampling each path and moving path samples locally toward locations with higher pressure, while ensuring that paths preserve the following properties: (i) they satisfy the maximum length constraint; (ii) to reach maximum coverage, they do not (self-)intersect; and (iii) they do not bend at sharp angles, as the design decisions of the algorithm stem from perceptual observations on smooth paths, and paths with sharp corners reach smaller coverage. To implement refinement, the goal and the properties are formulated as cost terms of an objective function, and a minimization algorithm is executed.

[00003]$P=\left\{x_i\right\}=\arg \underset{\left\{x_i\right\}}{\min }{c}_p+c_l+c_b$

[0067] Given paths with samples {x.sub.i ∈.sup.2}, we formulate a pressure intensity cost term as:

c.sub.p=−Σ.sub.ip*(x.sub.i)  (2)

[0068] The pressure intensity cost term (c.sub.p) is minimized as the path samples x.sub.i (i.e., points forming or defining the path) move to locations with higher pressure.

[0069] With a target path length L, and N samples per path, the target length is obtained if the length of each path segment is L/N. Then, a length cost term is formulated as:

[00004]$\begin{array}{cc}{c}_l={.Math.}_i{\left(.Math.x_{i+1}-x_i.Math.-\frac{L}{N}\right)}^2;& \left(3\right)\end{array}$

where x.sub.i and x.sub.i+1 are two consecutive path samples.

[0070] If two paths or two portions of a path get closer than the fall-off distance of focal points, σ, they stimulate overlapping skin areas. The result can be considered inefficient, as the covered skin area is larger if the paths move away. A (self-)intersection cost term that prevents path samples from getting too close is formulated as:

c.sub.i=Σ.sub.i,j max(σ−∥x.sub.j−x.sub.i∥, 0).sup.2;  (4)

where x.sub.iand x.sub.j are two non-consecutive path samples.

[0071] Finally, to favor low-curvature paths, a bending cost term is formulated as:

[00005]$\begin{array}{cc}{c}_b={.Math.}_i\mathrm{arc}{\tan }^2\frac{\left(x_{i+1}-x_i\right)\times \left(x_i-x_{i-1}\right)}{{\left(x_{i+1}-x_i\right)}^T\left(x_i-x_{i-1}\right)};& \left(5\right)\end{array}$

[0072] where x.sub.i−1, x.sub.i and x.sub.i+1 are three consecutive path samples.

[0073] In an embodiment, the paths are optimized by iterating steps of gradient descent of the four cost terms (although fewer cost terms or even additional cost terms may also be considered). For the pressure intensity term, a 2D grid is set with the pressure target values p*(x.sub.i), and bicubic interpolation is used to evaluate the pressure at sub-grid resolution and to compute robust gradients. Furthermore, a line search is applied to ensure that the step along the gradient reduces the cost. FIG. 5E shows an example of the computed paths 514 after refinement.

[0074] To account for the effective magnitude of the pressure field, a heuristic gain γ, as previously explained in equation (1), is preferably incorporated. The rendered pressure of a point on a path is set as p.sub.i=γ p*(x.sub.i), based on the target pressure field p*. FIG. 5F shows the reconstructed pressure field (result field 516) according to the quasi-static pressure model previously discussed in equation (1).

[0075] Once the 2D paths 514 are fully computed, they are lifted back to 3D for STM rendering on the ultrasound transducer array. This is achieved by undoing the projection of the hand voxels.

[0076] Each path is 140 mm long and is traversed in 20 ms at 7 m/s. The STM rendering API of an ultrasound transducer array updates every 1 ms a burst of X consecutive focal point positions (e.g., X being 40 in the Ultrahaptics STRATOS Explore USX). Therefore, each path is linearly up-sampled to 800 points spaced 0.175 mm. These points are fed in groups of 40 points to the STM rendering API of the ultrasound transducer array.

[0077] In accordance with a further aspect of the present invention there is provided a system for rendering tactile interaction using an ultrasound transducer array. The system comprises a data processing device configured to receive a target field (either a force field or a pressure field); compute one or more paths of pressure focal points of varying intensity such that the difference between a result field produced by the one or more paths of pressure focal points and the target field is minimized; and control an ultrasound transducer array in spatiotemporal modulation mode to simultaneously render the one or more computed one or more paths of pressure focal points (for instance, by computing spatiotemporal modulation commands corresponding to the 3D focal point paths and sending the commands to an ultrasound transducer array). The data processing device may be configured to implement the method as previously defined in any of the different embodiments.

[0078] FIG. 6 depicts a system 600 according to an exemplary embodiment. The system 600 comprises a data processing device 610, such as a computer with a processor and a GPU. The data processing device 610 manages the interaction of a user's hand 602 with gaseous fluid media in a virtual environment using an ultrasound transducer array 620. A tracking device 630 is configured to track the position of the body or a part of the body of the user. In the depicted example, the tracking device 630 is a hand-tracking device, such as a Leap Motion controller, that monitors the position of the user's hand 602.

[0079] Both the tracking device 630 and the ultrasound transducer array 620 are connected (e.g., via USB) to the computer running a simulation (for instance, a dynamic fluid simulation). The ultrasound transducer array 620 is actuated by a driver circuit 640 controlled by the computer. The system may comprise the driver circuit 640 configured to actuate the ultrasound transducer array 620. The system may also comprise the ultrasound transducer array 620. The system may comprise the tracking device 630 for tracking the position of the body or the part of the body of the user (in this particular case, for tracking the position of the user's hand 602).

[0080] The data processing device 610 implements the method as defined in any of the previous figures. The computed phase delays and amplitude values in STM mode are sent from the data processing device 610 to the driver circuit 640. The data processing device 610 may be configured to compute an interactive simulation (e.g., dynamic fluid simulation) and extract the target field 102 from the simulation. In this case, the target field 102 is generated from an interaction of a user's hand 602 with a virtual environment. The computer receives the position of the hand and updates the target field 102 generated on the user's hand 602 based on the hand position in the virtual environment.

[0081] The interactions produced during the virtual simulation may be shown on a display device 650, such as a monitor, configured to represent the interactive simulation in a virtual environment 704. FIG. 7 depicts the virtual representation (i.e., virtual hand 702) of the user's hand 602 interacting with virtual plumes of smoke 706 used in the dynamic fluid simulation. The plumes of smoke 706 generate a pressure target field 102 on the virtual hand 702 in the virtual environment 704.

[0082] The data processing device 610 computes the target field 102 and the paths 514 of pressure focal points of varying intensity that produce a best-matching pressure field on the user's hand, according to the method described. FIG. 8A depicts the pressure target field 102 produced on the surface of the hand by the columns of smoke of FIG. 7, which sets the target for the optimization algorithm of the method. FIG. 8B illustrates the reconstructed pressure field, result field 516, produced by STM rendering of dynamically optimized focal point paths.

[0083] In such interaction, haptic perception is dictated by a spatially and temporally varying pressure field on skin, which is used as target for the optimization algorithm. FIGS. 9 and 10 compare the reconstruction quality of the present rendering method vs. the rendering method of [8]. FIGS. 9A and 10A depict a 3D representation of a dynamic fluid simulation with one plume of smoke 706 and four plumes of smoke 706, respectively; FIGS. 9B and 10B depict the respective target field 102 generated on the virtual hand 702; FIGS. 9C and 10C show the respective result field 516 rendered by the present method whereas FIGS. 9D and 10D illustrates the reconstructed field 902 rendered by the AM-based method of [8]. These results show that the method of the present invention succeeds to provide larger and smoother coverage than the AM-based method of [8]. In particular, AM produces ambiguous results when rendering interaction with one wide plume (FIG. 9D) or with multiple thin plumes (FIG. 10D), while the present STM-based method does not suffer such ambiguity.