Calculating Normals analytically (Jacobian Matrix?)

OpenGL OpenGL: Advanced Coding
1

Hi I have this Gerstner function:

float PI = 3.14159265358979323846264;
P = pos;

         for (int i = 0; i < numWaves; ++i) { 
         float W = length(directions[i]);//wavelengths[i] ; //  wavelength 
         float frequency = 2.0 * PI / wavelengths[i]; // Frequency or k 

         float phase = speeds[i] * (2/W); // Phase 
         float Q = steepness[i]/(frequency * W * 4.0); // Steepness  / phase 





        P.x +=   (speeds[i] * (2/length))  * length *  directions[i].x * cos(frequency * dot(  directions[i].xy , pos.xy) +  (speeds[i] * (2/length)) * time); 
        P.y +=   (speeds[i] * (2/length))  * length *  directions[i].y * cos(frequency * dot(  directions[i].xy , pos.xy) +  (speeds[i] * (2/length)) * time); 
        P.z +=   length * cos(frequency * dot( directions[i].xy ,  pos.xy) + (speeds[i] * (2/length)) * time);

So I have a vector function f(x,y,z) -> (P.x,P.y,P.z).

I know I have to calculate the normal analytically as I’m displacing the vertex in the shader. I just have no idea how to go about it. I’ve tried getting the Jacobian
in Maple and Matlab, but aren’t getting the right results. Could someone with a little Calculus knowledge guide me in the right direction…?

Regards

2

You need to differentiate wrt pos.x and pos.y, then find their cross product.

The derivative of dot(directions[i].xy, pos.xy) wrt pos.x is frequency * directions[i].x, similarly wrt pos.y.

The derivative of cos(f(x)) wrt x is -sin(f(x))*f’(x).

So overall, you have:


Dx.x -= frequency * directions[i].x * (speeds[i] * (2/length)) * length * directions[i].x * sin(frequency * dot( directions[i].xy , pos.xy) + (speeds[i] * (2/length)) * time);
Dx.y -= frequency * directions[i].x * (speeds[i] * (2/length)) * length * directions[i].y * sin(frequency * dot( directions[i].xy , pos.xy) + (speeds[i] * (2/length)) * time);
Dx.z -= frequency * directions[i].x * length * sin(frequency * dot( directions[i].xy , pos.xy) + (speeds[i] * (2/length)) * time);

Dy.x -= frequency * directions[i].y * (speeds[i] * (2/length)) * length * directions[i].x * sin(frequency * dot( directions[i].xy , pos.xy) + (speeds[i] * (2/length)) * time);
Dy.y -= frequency * directions[i].y * (speeds[i] * (2/length)) * length * directions[i].y * sin(frequency * dot( directions[i].xy , pos.xy) + (speeds[i] * (2/length)) * time);
Dy.z -= frequency * directions[i].y * length * sin(frequency * dot( directions[i].xy , pos.xy) + (speeds[i] * (2/length)) * time);

...

normal = cross(Dx,Dy);

IOW, change cos() to sin(), multiply by frequency * directions[i].x (or .y), and negate (although that doesn’t actually matter because (-A)×(-B)=A×B).

Obviously, there’s a lot of potential to remove common subexpressions from that.

Also, for large numbers of waves, it’s common to use an inverse FFT rather than explicitly evaluating and summing sine waves. Sampling m waves at n points requires O(mn) time, while a FFT can do n waves at n points in O(nlog(n)).

3

BTW, I just noticed this, which may be a source of misunderstanding:

No you don’t; you have a function f(u,v) -> (P.x,P.y,P.z).

The difference doesn’t really matter much here; it just means that the Jacobian won’t be square (it will be 3x2). Although it may have caused some confusion if you were expecting the parameter (pos) to be 3D.