Perlin Noise: sharp contours between grid cells (Matlab)

Question

I believe this question is unrelated to this post.

I am trying to implement 2D classic Perlin noise in Matlab. I have been building off of this site as a reference. However, my program gives me "non-smooth" noise as you can see below. It appears that the local interpolation (inside each grid cell) works fine, but there is a continuity problem between grid cells. It appears to me that the noise is almost there, except for the sharpness, which I can't identify the cause of. If you look at the contour map, there is some matching between grid cells, so I am not (or I don't believe I am) picking the wrong grid points for the noise calculation.

.

Here is my code. I am not aiming for speed or efficiency at this point (it's very slow).

What is causing discontinuities along the grid?

function noise = getPerlinNoise(xPix, yPix, resolution, seed)

% allocate 2D noise map
noise = zeros(xPix, yPix);

% calculate grid spacing
gridSpacing = 1 / resolution;

% allocate gradient matrix
xGridSize = ceil(xPix/gridSpacing) + 1;
yGridSize = ceil(yPix/gridSpacing) + 1;
gradients = zeros(xGridSize, yGridSize, 2);

% seed PRNG for repeatable results
rng(seed);

% compute random gradient vectors (populate 'gradients' matrix)
for i = 1:xGridSize
    for j = 1:yGridSize

        % randomize gradients by choosing angle between 0 and 2*pi
        angle = rand * 2 * pi;

        % insert components into gradients matrix
        % Because point is on unit circle, no vector normalization is needed
        gradients(i, j, 1) = cos(angle);
        gradients(i, j, 2) = sin(angle);

    end
end

% actually calculate the noise
for i = 0:xPix - 1
    for j = 0:yPix - 1
        noise(i + 1, j + 1) = perlin(i, j);
    end
end

% ------
% END
% ------

%=====================================

    function val = lerp(a, b, w)
        % linearly interpolates between a and b with weight w
        val = a + w * (b - a);
    end

%=====================================

    function val = fade(t)
        % improved perlin fade function
        val = 6 * t^5 - 15 * t^4 + 10 * t^3;
    end

%=====================================


    function val = perlin(x, y)
        % computes perlin noise at pixel coordinates (x, y)

        % find grid points that a specific point will be located inside
        xi = floor(x / gridSpacing) + 1;     % +1 because Matlab arrays are 1-based 
        xf = xi + 1;
        yi = floor(y / gridSpacing) + 1;
        yf = yi + 1;

        % get coords of cell corners in pixel coordinates
        xCellMin = (xi - 1) * gridSpacing;
        xCellMax = xCellMin + gridSpacing;
        yCellMin = (yi - 1) * gridSpacing;
        yCellMax = yCellMin + gridSpacing;        

        % calculate distance vectors from grid points to current point
        A = [x - xCellMin, y - yCellMin];
        B = [x - xCellMax, y - yCellMin];
        C = [x - xCellMin, y - yCellMax];
        D = [x - xCellMax, y - yCellMax];        

        % calculate influence of gradient vectors by using dot product
        s = dot(A, [gradients(xi, yi, 1), gradients(xi, yi, 2)] );
        t = dot(B, [gradients(xi, yf, 1), gradients(xi, yf, 2)] );
        u = dot(C, [gradients(xf, yi, 1), gradients(xf, yi, 2)] );
        v = dot(D, [gradients(xf, yf, 1), gradients(xf, yf, 2)] );

        % interpolate to get final noise value at (x, y)
        sx = fade((x - xCellMin) / gridSpacing);
        sy = fade((y - yCellMin) / gridSpacing);

        a = lerp(s, t, sx);
        b = lerp(u, v, sx);
        val = lerp(a, b, sy);      

    end
end

==============================================

Corrected Code

(thanks to MarkV)

    % calculate influence of gradient vectors by using dot product
    s = dot(A, [gradients(xi, yi, 1), gradients(xi, yi, 2)] );
    t = dot(B, [gradients(xi, yf, 1), gradients(xi, yf, 2)] );
    u = dot(C, [gradients(xf, yi, 1), gradients(xf, yi, 2)] );
    v = dot(D, [gradients(xf, yf, 1), gradients(xf, yf, 2)] );

should be changed to

    % calculate influence of gradient vectors by using dot product
    s = dot(A, [gradients(xi, yi, 1), gradients(xi, yi, 2)] );
    t = dot(B, [gradients(xf, yi, 1), gradients(xf, yi, 2)] );
    u = dot(C, [gradients(xi, yf, 1), gradients(xi, yf, 2)] );
    v = dot(D, [gradients(xf, yf, 1), gradients(xf, yf, 2)] );

Answer 1

It seems to me like you switched B & C here.

% calculate influence of gradient vectors by using dot product
s = dot(A, [gradients(xi, yi, 1), gradients(xi, yi, 2)] );
t = dot(C, [gradients(xi, yf, 1), gradients(xi, yf, 2)] ); % use C
u = dot(B, [gradients(xf, yi, 1), gradients(xf, yi, 2)] ); % use B
v = dot(D, [gradients(xf, yf, 1), gradients(xf, yf, 2)] );