2D Water Shader

This shader distorts textures behind it according to a noise map and also creates a nice blue/green color gradient. Requires a noise map but not a base texture.

// Upgrade NOTE: replaced 'mul(UNITY_MATRIX_MVP,*)' with 'UnityObjectToClipPos(*)'

 Shader "Custom/WaterShader2D" 
 {
     Properties 
     {
         _MainTex ("Base (RGB) Trans (A)", 2D) = "white" {}
         _Colour ("Colour", Color) = (1,1,1,1)
 
         _BumpMap ("Noise text", 2D) = "bump" {}
         _Magnitude ("Magnitude", Range(0,1)) = 0.5
		 _Blue ("Blue", Range(0,7)) = 2
		 _Green ("Green", Range(0,7)) = 2
		 _MaxAlpha("Alpha Clamp", Range(0,1)) = 0.5
     }
     
     SubShader
     {
         Tags {"Queue"="Transparent" "IgnoreProjector"="True" "RenderType"="Opaque"}
         ZWrite On Lighting Off Cull Off Fog { Mode Off } Blend One Zero
 
         GrabPass { "_GrabTexture" }
         
         Pass 
         {
             CGPROGRAM
             #pragma vertex vert
             #pragma fragment frag
             #include "UnityCG.cginc"
 
             sampler2D _GrabTexture;
 
             sampler2D _MainTex;
             fixed4 _Colour;
 
             sampler2D _BumpMap;
             float  _Magnitude;
			 float _Green;
			 float _Blue;
			 float _MaxAlpha;
			 
 
             struct vin_vct
             {
                 float4 vertex : POSITION;
                 float4 color : COLOR;
                 float2 texcoord : TEXCOORD0;
				 //float2 uv : TEXCOORD0;
             };
 
             struct v2f_vct
             {
                 float4 vertex : POSITION;
                 fixed4 color : COLOR;
                 float2 texcoord : TEXCOORD0;
				 //float2 uv : TEXCOORD0;
                 float4 uvgrab : TEXCOORD1;
             };
 
             // Vertex function 
             v2f_vct vert (vin_vct v)
             {
                 v2f_vct o;
                 o.vertex = UnityObjectToClipPos(v.vertex);
                 o.color = v.color;
 
                 o.texcoord = v.texcoord;
				 //o.uv = TRANSFORM_TEX(v.uv, _MainTex);
                 o.uvgrab = ComputeGrabScreenPos(o.vertex);
                 return o;
             }
 
             // Fragment function
             half4 frag (v2f_vct i) : COLOR
             {
                 half4 mainColour;
                 mainColour.a = 1;
				 mainColour.b = 1*_Blue*i.texcoord.y;
				 mainColour.g = 1*_Green*i.texcoord.y*i.texcoord.y;
				 mainColour.r = 0;
                 half4 bump = tex2D(_BumpMap, i.texcoord);
                 half2 distortion = UnpackNormal(bump).rg;
 
                 i.uvgrab.xy += distortion * _Magnitude*.02;
 
                 fixed4 col = tex2Dproj( _GrabTexture, UNITY_PROJ_COORD(i.uvgrab));
				 col.a = 0;
                 return col * mainColour ;
             }
         
             ENDCG
         } 
     }
 }

 

Modeling a Boat’s Wake – Part 4 – buoyant sphere interaction

 

After a pretty substantial amount of debugging (not to mention pages of math) – I have at least some working version of the interactive wave algorithm along with a super basic buoyancy model.

Ta Da! At present, there is still much work to do in terms of optimization and tweaking of wave parameters.  The dampening, gravity and time step parameters from Part 2 are extremely sensitive.  Also, the convolution applied from Tessendorf’s paper is written asymmetrically for some reason in k and l.  This resulted in waves not propagating quite radially.  When I made the convolution symmetric the computation become much more expensive. I need to do some more digging.

Modeling a Boat’s Wake – Part 3 – Code

Finally, we can write some code and put the math in part 2 to use. The first thing we need to do is write a few math functions.  First, lets take care of that Bessel function. For this, I pulled out my good old copy of Numerical Recipes and flipped to their special functions chapter.  For legal reasons I’m not going to post the code as its more or less verbatim taken from the NR pages.

    float bessj0(float x)
    //Returns the Bessel function J0(x) for any real x

 

I will say it was necessary to substitute NR’s dummy math functions with unity’s Mathf functions. And as such, I had to cast a few of the doubles used by the NR method as floats.  Though the bessj0 takes and returns a float, it uses doubles for most of its calculations. With a working Bessel function in my pocket, I can now start to write some methods to solve for the convolution kernel in Part 2.

    float G0()
    {
        float ans = 0;
        
        float qnsq;
        for (int i = 1; i<n;i ++)
        {
            qnsq = i*i*dq*dq;;
            ans = ans + qnsq * Mathf.Exp(-sigma * qnsq);
        }
        return ans;
    }

 

The above simple function generates the normalization constant, while :

    float[,] kernel()
    {
        float[,] ans = new float[P + 1, P + 1];
        Vector2 r;
        float qnsq;
        float tempAns;
        float Gnorm = G0();
       
        for (int k=0; k<P+1; k++)
        {
            for (int l = 0; k < P + 1; l++)
            {
                tempAns = 0;
                for (int i =1; i<n; i++)
                {
                    r.x = k;
                    r.y = l;
                    qnsq = Mathf.Pow(i * dq, 2);
                    tempAns = tempAns + qnsq * Mathf.Exp(-sigma * qnsq) * bessj0(i * dq * r.magnitude);
                }
                ans[k, l] = tempAns / Gnorm;

            }
        }
        return ans;

    }

 

This kernel function should calculate G(k,l).  Next we need to perform the ‘vertical derivative’.

    void convolveVerticalDerivative()
    {
        

        for (int i = 0;i< xRes;i++)
        {
            for (int j = 0; j < yRes; j++)
            {
                float temp = 0;
                for (int k = 0; k < P+1; k++)
                {
                    for (int l = k+1; l < P + 1; l++)
                    {
                        temp = temp + kernel[k, l] * (heightMap[i + k, j + l] + heightMap[i - k, j - l] + heightMap[i + k, j - l] + heightMap[i - k, j + l]);
                    }

                }
                verticalDerivative[i, j] = heightMap[i, j] + temp;

            }

        }

    }

 

The quad-nested for loop has me a bit concerned performance wise, especially for something that is called every frame.  However, if P is small and the heightmap grid isn’t outrageously large, this should be manageable. If not we may have to be smarter. I should mention that I’m doing something incredibly stupid, which is writing all of this code without doing intermediate testing – I am sure to pay for this, or maybe it’ll all just work?  We’ll see if my laziness pays off.

Now finally, the cou de gras, the actual beef of the code – we will write something to solve the remainder of the PDE and update the height map.

    void updateHeightMap()
    {
        updateSource();
        convolveVerticalDerivative();
        float temp;
        for (int i = 0; i < xRes; i++)
        {
            for (int j = 0; j < yRes; j++)
            {
                
                temp = heightMap[i, j];
                heightMap[i, j] = heightMap[i, j] * (2.0f - alpha * dt) / (1.0f + alpha * dt) - prevHeightMap[i, j] / (1.0f + alpha * dt) - verticalDerivative[i, j] * g * dt * dt / (1.0f + alpha * dt);
                prevHeightMap[i, j] = temp;
            }
        }

    }

 

So, theoretically, our interactive water code should be working. The small parts seem to be working in the limited testing I’ve done, but for anything real to happen I need to interface this code to  Unity mesh.  And for it to be ‘interactive’, I need some way to edit the heightMap to create disturbances.  Also, I need to call these functions in the standard Unity start and update functions

    void Start()
    {
        calculateKernel();
    }

    void FixedUpdate()
    {
        updateHeightMap();

    }

    void Update()
    {
        water.setHeightmap(heightMap);
    }

    void updateSource()
    {
        for (int i = 0; i < xRes; i++)
        {
            for (int j = 0; j < yRes; j++)
            {
                heightMap[i, j] = heightMap[i, j] + source[i, j];
            }
        }
    }

 

In the start function we simply initialize the convolution kernel, in fixed update we advance the height map – and in update we will update our water mesh with the new heightMap.  In unity, FixedUpdate is called at a fixed interval in game time. For this reason, we should do all of our physics in FixedUpdate.  Update however, is a good place to do mesh and other graphics updates.  Finally, we have updateSource – this function is not a Unity function but rather a function which we will call in updateHeightMap which will allow us to interact with the water.  For example, we can have the source array represent water displaced momentarily by an object.

Well that’s it for the bulk of the actual code – next time we’ll work on testing and full integration with Unity.

 

 

Modeling a Boat’s Wake – Part 2 – Maths

Equations of Motion

In order to model hydrodynamic surface waves, I am following a general model for interactive water graphics based on Bernoulli’s equation. This approach, along with the basic equations I will use, are those outlined in Tessendorf’s paper ‘Interactive Water Surfaces’.[1]  The general algorithm he develops is what I will implement.  I will try to supplement his derivation rather than directly list his steps.  The general equation of motion used is a linearized Bernoulli’s equation as is below. [1]

The first two terms will be solved using rather straight forward finite difference methods.  The Laplacian term on the RHS will be solved in a similar, yet more complex manner: use of a convolution kernel.

Convolution Kernel as a Linear Operator

Written more explicitly, the calculus of interest is as follows

This term, with the strange negative root of a Laplacian operator, apparently accounts for conservation of mass and a gravitational restoring force.  The g is a gravitational constant while h is our 2D time dependent height map. I would be interested to see how this strange operator appears but we’ve got more pressing issues.  To numerically calculate this quantity efficiently we are going to need to make use of a convolution kernel – and if you don’t do image processing on a regular basis like me you may need a bit of background to swallow how they work.   The following explanation of convolution follows very closely this lecture/paper on image processing [2] .

So lets start with some hopefully familiar ground, good ole’ standard definition of a partial derivative:

Now numerically we could do this several ways, most common are simple two point methods and higher order methods using several weighted points.  We will see that one advantage of using convolution kernels is that changing  this method numerically is simple and the general code is very generic.  Anyhow, numerically we could estimate this partial derivative roughly as:

Where i and j are integer indices corresponding to x = i*dx, y = j*dy and dx =1, AKA finite difference. Similarly we could do a centered finite difference

Or we could get even fancier a do a five point method

Now we can begin to see that this could be rewritten as a second function modifying or convolving our function F.  This convolution function can as Tessendorf points out, perform any linear operator on F.  But lets not get too far ahead, lets rewrite this five point stencil method as

or

Where G(k) is simply a function defining the coefficients above:  [1/12, -8/12, 0, 8/12, -1/12]. In this context, the space over which the convolution function operates (from k=-2 to 2 in this case), is known as the kernel.

Okay now, so what about that somewhat comparatively nasty square root of negative 2D Laplacian? Well Tessendorf of course applies an appropriate convolution kernel.  How he chose this convolution function is currently a mystery to me, I suppose similar to how the exact coefficients in the five point stencil method were chosen – which I am also in the dark on.  If anyone has some insight, do let me know.  For the mean time, I will assume that these weighting functions are probably effective depending on the function you are differentiating – and that Tessendorf chose a good function for this application.  Again, if anyone knows anything about this decision I am truly curious. Anyhow, here we go. Tessendorf writes the operator acting on h(x,y,t) as

where

Jo is a Bessel function of the first kind, N is some integer such that N>> n , sigma is chosen so that the series is convergent, Tessendorf uses 1,  , and r is simply

The normalization constant, Go, is given by

With this ‘convolution kernel’ now known and somewhat demystified (at least for me), we can move on to solving the full differential equation.  There is a bit of a numerical trick in the paper in which Tessendorf exploits the symmetry of the convolution kernel to drastically reduce the number of operations.

Finite Difference – Acceleration and Velocity terms

Now on the left hand side of the linearized Bernoulli’s equation we have our vertical acceleration and a dampening term, both derivatives in time.  The terms are as follows

for the dampening term on the right, we can simply write it as

While the second derivative on the left is a bit more complicated.  Following [1] we will use a second order central difference.

All together our numerical differential equation looks like

In order to advance our solution to the next time – or rather in order to propagate the wave, we need to solve for 

Someone who is really paying attention (I would be super surprised), would notice that my final equation above has a slight difference than in [1].  I have 2+dt alpha, while [1] has 2-dt alpha in the first term on the RHS.  I may be wrong, but I double and triple checked this and got the same result each time.  It may be that I found a mistake – seems unlikely, I would follow [1] if I were you and use 2-dt alpha, but I have to see this thing through and try +.

So that’s it, in terms of Math, now we just need to write some code!

[1] Jerry Tessendorf, Interactive Water Surfaces

[2] Tomasi, Carlos, Convolution, Smoothing, and Image Derivatives , https://www.cs.duke.edu/courses/spring03/cps296.1/handouts/Image%20Processing.pdf

 

Modeling a Boat’s Wake – Part 1 – Background

The Kelvin Wake

The goal of my first ‘big’ Unity3D project is to create some realistic sailboat physics including some relatively accurate hydrodynamics and aerodynamics.  Currently, one of my focuses is on realistic water graphics/physics.  One effect I absolutely want to be able to model the Kelvin Wake, which I’m sure looks familiar ( picture from the good folks over at http://www.boatdesign.net/forums/)

The general equations necessary to accurately model a Kelvin Wake turn out to be relatively difficult to solve – at least with any computational efficiency.  My first idea to cheaply recreate this effect using a Unity mesh was to import an image similar to the above using a stand alone script which could be called as a function.  I would have to convert the image to a scaled 2D matrix varying from 0-1 or similar.  This process would not be unlike what may people do to import height maps for Terrains.  Then I could modify the function based on my objects current and historical velocity vector.  This all sounds pretty do-able and maybe even efficient if the effect could be created in a believable or even good looking fashion.  That is until I saw this, a working demo of interactive water which could potentially recreate the Kelvin Wake along with a multitude of other realistic effects.

After a bit more research I found, not surprisingly, that people have written a great deal on interactive water computer graphics for precisely the applications I am seeking. Below is a brief list of what I’ve been reading :

[1] J. Loviscach, A Convolution-Based Algorithm for Animated Water Waves, EUROGRAPHICS 2002

[2] Jerry Tessendorf. Interactive water surfaces. In Game Programming Gems 4. Charles River Media, 2004

[3] Jerry Tessendorf. Simulating ocean surface. SIGGRAPH 2004 Course Notes, 2004.

http://tessendorf.org/reports.html

This Tessendorf fellow seems to be quite prolific realm of simulating water dynamics for game programming.  The second reference there seems to have the simplest algorithm (iWave) to follow for basic interactive water surfaces.   – I’m gonna do a bit more reading and head to bed.  Cheers.

Hello Word

This place, as it stands, is meant to serve as a working ‘lab notebook’ for my adventures in learning Unity3D along with the basics of game development, computer graphics and c#.  My background is in physics so I have some ground to stand on and have written plenty of code in my life.  Game development isn’t too far a cry from computational physics, but it sure seems like more fun. Anyhow, take everything here with a grain of salt and thanks for stopping by.