Many years ago (mid-80s), I was working at a company that had a Silicon Graphics workstation.  Among a handful of demos designed to show off the SGI machine’s high-end graphics was a simulation of wave propagation in a little wireframe mesh.  It was great fun to play with by changing the height of points in the mesh and then letting the simulation run.  And the SGI machine was fast enough that the resulting animation was just mesmerizing.

Recreating this water simulation in WPF seemed like a nice way to learn a little more about 3D graphics in WPF.  (The end result is here).

The first step was to find an algorithm that simulates wave propagation through water.  It turns out that there is a very simple algorithm that achieves the desired effect simply by taking the average height of neighboring points.  The basic algorithm is described in detail in this article on 2D Water.  The same algorithm is also described in The Water Effect Explained.

The next step is to set up the 3D viewport and its constituent elements.  I used two different directional lights, to create more contrast on the surface of the water, as well as defining both diffuse and specular material properties for the surface of the water.

Here is the relevant XAML.  Note that meshMain is the mesh that will contain the surface of the water.

        <Viewport3D Name="viewport3D1" Margin="0,8.181,0,0" Grid.Row="1">
            <Viewport3D.Camera>
                <PerspectiveCamera x:Name="camMain" Position="48 7.8 41" LookDirection="-48 -7.8 -41" FarPlaneDistance="100" UpDirection="0,1,0" NearPlaneDistance="1" FieldOfView="70">

                </PerspectiveCamera>
            </Viewport3D.Camera>
            <ModelVisual3D x:Name="vis3DLighting">
                <ModelVisual3D.Content>
                    <DirectionalLight x:Name="dirLightMain" Direction="2, -2, 0"/>
                </ModelVisual3D.Content>
            </ModelVisual3D>
            <ModelVisual3D>
                <ModelVisual3D.Content>
                    <DirectionalLight Direction="0, -2, 2"/>
                </ModelVisual3D.Content>
            </ModelVisual3D>
            <ModelVisual3D>
                <ModelVisual3D.Content>
                    <GeometryModel3D x:Name="gmodMain">
                        <GeometryModel3D.Geometry>
                            <MeshGeometry3D x:Name="meshMain" >
                            </MeshGeometry3D>
                        </GeometryModel3D.Geometry>
                        <GeometryModel3D.Material>
                            <MaterialGroup>
                                <DiffuseMaterial x:Name="matDiffuseMain">
                                    <DiffuseMaterial.Brush>
                                        <SolidColorBrush Color="DarkBlue"/>
                                    </DiffuseMaterial.Brush>
                                </DiffuseMaterial>
                                <SpecularMaterial SpecularPower="24">
                                    <SpecularMaterial.Brush>
                                        <SolidColorBrush Color="LightBlue"/>
                                    </SpecularMaterial.Brush>
                                </SpecularMaterial>
                            </MaterialGroup>
                        </GeometryModel3D.Material>
                    </GeometryModel3D>
                </ModelVisual3D.Content>
            </ModelVisual3D>
        </Viewport3D>

Next, we create a WaveGrid class that implements the basic algorithm described above.  The basic idea is that we maintain two separate buffers of mesh data—one representing the current state of the water and one the prior state.  WaveGrid stores this data in two Point3DCollection objects.  As we run the simulation, we alternate which buffer we’re writing into and attach our MeshGeometry3D.Positions property to the most recent buffer.  Note that we’re only changing the vertical height of the points—which is the Y value.

WaveGrid also builds up the triangle indices for the mesh, in an Int32Collection which will also get connected to our MeshGeometry3D.

All of the interesting stuff happens in ProcessWater.  This is where we implement the smoothing algorithm described in the articles.  Since I wanted to fully animate every point in the mesh, I processed not just the internal points that have four neighboring points, but the points along the edge of the mesh, as well.  As we add in height values of neighboring points, we keep track of how many neighbors we found, so that we can do the averaging properly.

The final value for each point is a function of both the smoothing (average height of your neighbors) and the “velocity”, which is basically—how far from equilibrium was the point during the last iteration?  We also then apply a damping factor, since waves will gradually lose their amplitude.

Here’s the complete code for the WaveGrid class:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows.Media;
using System.Windows.Media.Media3D;

namespace WaveSim
{
    class WaveGrid
    {
        // Constants
        const int MinDimension = 5;    
        const double Damping = 0.96;
        const double SmoothingFactor = 2.0;     // Gives more weight to smoothing than to velocity

        // Private member data
        private Point3DCollection _ptBuffer1;
        private Point3DCollection _ptBuffer2;
        private Int32Collection _triangleIndices;

        private int _dimension;

        // Pointers to which buffers contain:
        //    - Current: Most recent data
        //    - Old: Earlier data
        // These two pointers will swap, pointing to ptBuffer1/ptBuffer2 as we cycle the buffers
        private Point3DCollection _currBuffer;
        private Point3DCollection _oldBuffer;

        /// <summary>
        /// Construct new grid of a given dimension
        /// </summary>
        ///
<param name="Dimension"></param>
        public WaveGrid(int Dimension)
        {
            if (Dimension < MinDimension)
                throw new ApplicationException(string.Format("Dimension must be at least {0}", MinDimension.ToString()));

            _ptBuffer1 = new Point3DCollection(Dimension * Dimension);
            _ptBuffer2 = new Point3DCollection(Dimension * Dimension);
            _triangleIndices = new Int32Collection((Dimension - 1) * (Dimension - 1) * 2);

            _dimension = Dimension;

            InitializePointsAndTriangles();

            _currBuffer = _ptBuffer2;
            _oldBuffer = _ptBuffer1;
        }

        /// <summary>
        /// Access to underlying grid data
        /// </summary>
        public Point3DCollection Points
        {
            get { return _currBuffer; }
        }

        /// <summary>
        /// Access to underlying triangle index collection
        /// </summary>
        public Int32Collection TriangleIndices
        {
            get { return _triangleIndices; }
        }

        /// <summary>
        /// Dimension of grid--same dimension for both X & Y
        /// </summary>
        public int Dimension
        {
            get { return _dimension; }
        }

        /// <summary>
        /// Set center of grid to some peak value (high point).  Leave
        /// rest of grid alone.  Note: If dimension is even, we're not
        /// exactly at the center of the grid--no biggie.
        /// </summary>
        ///
<param name="PeakValue"></param>
        public void SetCenterPeak(double PeakValue)
        {
            int nCenter = (int)_dimension / 2;

            // Change data in oldest buffer, then make newest buffer
            // become oldest by swapping
            Point3D pt = _oldBuffer[(nCenter * _dimension) + nCenter];
            pt.Y = (int)PeakValue;
            _oldBuffer[(nCenter * _dimension) + nCenter] = pt;

            SwapBuffers();
        }

        /// <summary>
        /// Leave buffers in place, but change notation of which one is most recent
        /// </summary>
        private void SwapBuffers()
        {
            Point3DCollection temp = _currBuffer;
            _currBuffer = _oldBuffer;
            _oldBuffer = temp;
        }

        /// <summary>
        /// Clear out points/triangles and regenerates
        /// </summary>
        ///
<param name="grid"></param>
        private void InitializePointsAndTriangles()
        {
            _ptBuffer1.Clear();
            _ptBuffer2.Clear();
            _triangleIndices.Clear();

            int nCurrIndex = 0;     // March through 1-D arrays

            for (int row = 0; row < _dimension; row++)
            {
                for (int col = 0; col < _dimension; col++)
                {
                    // In grid, X/Y values are just row/col numbers
                    _ptBuffer1.Add(new Point3D(col, 0.0, row));

                    // Completing new square, add 2 triangles
                    if ((row > 0) && (col > 0))
                    {
                        // Triangle 1
                        _triangleIndices.Add(nCurrIndex - _dimension - 1);
                        _triangleIndices.Add(nCurrIndex);
                        _triangleIndices.Add(nCurrIndex - _dimension);

                        // Triangle 2
                        _triangleIndices.Add(nCurrIndex - _dimension - 1);
                        _triangleIndices.Add(nCurrIndex - 1);
                        _triangleIndices.Add(nCurrIndex);
                    }

                    nCurrIndex++;
                }
            }

            // 2nd buffer exists only to have 2nd set of Z values
            _ptBuffer2 = _ptBuffer1.Clone();
        }

        /// <summary>
        /// Determine next state of entire grid, based on previous two states.
        /// This will have the effect of propagating ripples outward.
        /// </summary>
        public void ProcessWater()
        {
            // Note that we write into old buffer, which will then become our
            //    "current" buffer, and current will become old. 
            // I.e. What starts out in _currBuffer shifts into _oldBuffer and we
            // write new data into _currBuffer.  But because we just swap pointers,
            // we don't have to actually move data around.

            // When calculating data, we don't generate data for the cells around
            // the edge of the grid, because data smoothing looks at all adjacent
            // cells.  So instead of running [0,n-1], we run [1,n-2].

            double velocity;    // Rate of change from old to current
            double smoothed;    // Smoothed by adjacent cells
            double newHeight;
            int neighbors;

            int nPtIndex = 0;   // Index that marches through 1-D point array

            // Remember that Y value is the height (the value that we're animating)
            for (int row = 0; row < _dimension ; row++)
            {
                for (int col = 0; col < _dimension; col++)
                {
                    velocity = -1.0 * _oldBuffer&#91;nPtIndex&#93;.Y;     // row, col
                    smoothed = 0.0;

                    neighbors = 0;
                    if (row > 0)    // row-1, col
                    {
                        smoothed += _currBuffer[nPtIndex - _dimension].Y;
                        neighbors++;
                    }

                    if (row < (_dimension - 1))   // row+1, col
                    {
                        smoothed += _currBuffer&#91;nPtIndex + _dimension&#93;.Y;
                        neighbors++;
                    }

                    if (col > 0)          // row, col-1
                    {
                        smoothed += _currBuffer[nPtIndex - 1].Y;
                        neighbors++;
                    }

                    if (col < (_dimension - 1))   // row, col+1
                    {
                        smoothed += _currBuffer&#91;nPtIndex + 1&#93;.Y;
                        neighbors++;
                    }

                    // Will always have at least 2 neighbors
                    smoothed /= (double)neighbors;

                    // New height is combination of smoothing and velocity
                    newHeight = smoothed * SmoothingFactor + velocity;

                    // Damping
                    newHeight = newHeight * Damping;

                    // We write new data to old buffer
                    Point3D pt = _oldBuffer&#91;nPtIndex&#93;;
                    pt.Y = newHeight;   // row, col
                    _oldBuffer&#91;nPtIndex&#93; = pt;

                    nPtIndex++;
                }
            }

            SwapBuffers();
        }
    }
}
&#91;/sourcecode&#93;

Finally, we need to hook everything up.  When our main window fires up, we create an instance of <strong>WaveGrid </strong>and set the center point in the grid to some peak value.  When we start the animation, this higher point will fall and trigger the waves.

We do all of the animation in the <strong>CompositionTarget.Rendering </strong>event handler.  This is the recommended spot to do custom animations in WPF, as opposed to doing the animation in some timer Tick event.  (<em>Windows Presentation Foundation Unleashed</em>, Nathan, pg 470).

When you attach a handler to the <strong>Rendering </strong>event, WPF just continues rendering frames indefinitely.  One problem is that the handler will get called for every frame rendered, which turns out to be too fast for our water animation.  To get the water to look right, we keep track of the time that we last rendered a frame and then wait a specified number of milliseconds before rendering another.

Here is the full source code for Window1.xaml.cs:



using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Data;
using System.Windows.Documents;
using System.Windows.Input;
using System.Windows.Media;
using System.Windows.Media.Media3D;
using System.Windows.Media.Imaging;
using System.Windows.Navigation;
using System.Windows.Shapes;
using System.Windows.Threading;

namespace WaveSim
{
    /// <summary>
    /// Interaction logic for Window1.xaml
    /// </summary>
    public partial class Window1 : Window
    {
        private Vector3D zoomDelta;

        private WaveGrid _grid;
        private bool _rendering;
        private double _lastTimeRendered;
        private double _firstPeak = 6.5;

        // Values to try:
        //   GridSize=20, RenderPeriod=125
        //   GridSize=50, RenderPeriod=50
        private const int GridSize = 50;   
        private const double RenderPeriodInMS = 50;    

        public Window1()
        {
            InitializeComponent();

            _grid = new WaveGrid(GridSize);        // 10x10 grid
            slidPeakHeight.Value = _firstPeak;
            _grid.SetCenterPeak(_firstPeak);
            meshMain.Positions = _grid.Points;
            meshMain.TriangleIndices = _grid.TriangleIndices;

            // On each WheelMouse change, we zoom in/out a particular % of the original distance
            const double ZoomPctEachWheelChange = 0.02;
            zoomDelta = Vector3D.Multiply(ZoomPctEachWheelChange, camMain.LookDirection);
        }

        private void Window_MouseWheel(object sender, MouseWheelEventArgs e)
        {
            if (e.Delta > 0)
                // Zoom in
                camMain.Position = Point3D.Add(camMain.Position, zoomDelta);
            else
                // Zoom out
                camMain.Position = Point3D.Subtract(camMain.Position, zoomDelta);
            Trace.WriteLine(camMain.Position.ToString());
        }

        // Start/stop animation
        private void btnStart_Click(object sender, RoutedEventArgs e)
        {
            if (!_rendering)
            {
                _grid = new WaveGrid(GridSize);        // 10x10 grid
                _grid.SetCenterPeak(_firstPeak);
                meshMain.Positions = _grid.Points;

                _lastTimeRendered = 0.0;
                CompositionTarget.Rendering += new EventHandler(CompositionTarget_Rendering);
                btnStart.Content = "Stop";
                slidPeakHeight.IsEnabled = false;
                _rendering = true;
            }
            else
            {
                CompositionTarget.Rendering -= new EventHandler(CompositionTarget_Rendering);
                btnStart.Content = "Start";
                slidPeakHeight.IsEnabled = true;
                _rendering = false;
            }
        }

        void CompositionTarget_Rendering(object sender, EventArgs e)
        {
            RenderingEventArgs rargs = (RenderingEventArgs)e;
            if ((rargs.RenderingTime.TotalMilliseconds - _lastTimeRendered) > RenderPeriodInMS)
            {
                // Unhook Positions collection from our mesh, for performance
                // (see http://blogs.msdn.com/timothyc/archive/2006/08/31/734308.aspx)
                meshMain.Positions = null;

                // Do the next iteration on the water grid, propagating waves
                _grid.ProcessWater();

                // Then update our mesh to use new Z values
                meshMain.Positions = _grid.Points;

                _lastTimeRendered = rargs.RenderingTime.TotalMilliseconds;
            }
        }

        private void slidPeakHeight_ValueChanged(object sender, RoutedPropertyChangedEventArgs<double> e)
        {
            _firstPeak = slidPeakHeight.Value;
            _grid.SetCenterPeak(_firstPeak);
        }
    }
}

The end result is pretty satisfying—a nice smooth animation of a series of ripples propagating out from the initial disturbance.  You can install and run the simulation by clicking here.  Note that you can zoom in/out using the mouse wheel.

We could extend this example in several different ways:

• Render the surface of the water in a more lifelike way—e.g. glassy, with reflections.
• Add simple controls to change the viewpoint or to rotate the mesh itself
• Add knobs for playing with things like Damping and SmoothingFactor
• Add ability to “grab” points in the mesh with the mouse and manually move them up/down
• Raindrop simulation—just add timer that introduces new random peaks, representing raindrops
• Antialiasing–also consider diagonally adjacent points as neighbors, but adjust by weighting factor when averaging