PDA

View Full Version : Calculating Normals analytically (Jacobian Matrix?)

Moetman
07-02-2017, 06:36 PM
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

GClements
07-03-2017, 06:41 AM
Could someone with a little Calculus knowledge guide me in the right direction..?
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(m*n) time, while a FFT can do n waves at n points in O(n*log(n)).

GClements
07-03-2017, 06:49 AM
BTW, I just noticed this, which may be a source of misunderstanding:

So I have a vector function f(x,y,z) -> (P.x,P.y,P.z).
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.