"Blocky" Perlin noise

Question

I've recently been trying to implement a Perlin Noise generator in C (based off Ken Perlin's website, using the SDL library as a screen output), but the output shows that the edges between interpolation blocks are not continuous or smooth - the interpolation blocks do really manifest as blocks.

I've tried four kinds of interpolations, and all the "smooth" ones look about the same; only cosine looks (very) slightly better and straight linear looks horrible by comparison. (below are Cosine and Linear)

Ironically, if making a fractal sum of noises (my ultimate purpose of this), linear blows away the smooth interpolations in terms of "blockyness" and actually looks almost fine.

I'm pretty sure that there's something I'm missing in my code or doing it wrong, but I can't seem to find it.

Any suggestions of what (or what conditions) may be causing these block artifacts?

For reference, my current code below:

#include<stdio.h>
#include<math.h>
#include<SDL/SDL.h>

void normalize3(float *vec3){
    float distX=0,distY=0,distZ=0;
    distX=vec3[0];
    distX*=distX;
    distY=vec3[1];
    distY*=distY;
    distZ=vec3[2];
    distZ*=distZ;
    float dist=sqrtf(distX+distY+distZ);
    vec3[0]/=dist;
    vec3[1]/=dist;
    vec3[2]/=dist;
}

float sinterpolate(float scale){
    //return scale*scale*(3.0-2*scale); //Classic 3*t^2-2*t^3

    /*float t=scale*scale;
    float u=t*t;
    return (6.0*u*scale-15.0*u+10.0*t*scale);*/ //Improved 6*t^5-15*t^4+10*t^3

    return (0.5-cosf(scale*M_PI)/2.0); //Straight cosine interpolation
}

float linterpolate(float a,float b,float scale){
    return a+scale*(b-a);
}

float noise3(float *vec3,float *grads,Uint8 *perms){
    vec3[0]=fmodf(vec3[0],256.0);
    vec3[1]=fmodf(vec3[1],256.0);
    vec3[2]=fmodf(vec3[2],256.0);
    Uint8 ivec3[3];

    float relPos[3],temp;
    float cube[2][2][2];
    Uint8 index;

    //One loop for each dimension of noise.
    for(int x=0;x<2;x++){
        ivec3[0]=vec3[0];
        ivec3[0]+=x;
        relPos[0]=vec3[0]-ivec3[0];
        for(int y=0;y<2;y++){
            ivec3[1]=vec3[1];
            ivec3[1]+=y;
            relPos[1]=vec3[1]-ivec3[1];
            for(int z=0;z<2;z++){
                ivec3[2]=vec3[2];
                ivec3[2]+=z;
                relPos[2]=vec3[2]-ivec3[2];

                index=ivec3[0]+perms[ivec3[1]+perms[ivec3[2]]];

                temp=relPos[0]*grads[3*index];
                temp+=relPos[1]*grads[3*index+1];
                temp+=relPos[2]*grads[3*index+2]; //The gradient's dot product
                                                  //with respect to the point
                                                  //being analyzed

                cube[x][y][z]=temp;
            }
        }
    }

    ivec3[0]--;
    ivec3[1]--;
    ivec3[2]--;
    relPos[0]=vec3[0]-ivec3[0];
    relPos[1]=vec3[1]-ivec3[1];
    relPos[2]=vec3[2]-ivec3[2];
    relPos[0]=sinterpolate(relPos[0]);  //Comment these
    relPos[1]=sinterpolate(relPos[1]);  //if you want
    relPos[2]=sinterpolate(relPos[2]);  //Linear Interpolation.


    return linterpolate(linterpolate(linterpolate(cube[0][0][0],cube[0][0][1],relPos[2]),linterpolate(cube[0][8][0], cube[0][9][1],relPos[2]),relPos[1]),linterpolate(linterpolate(cube[1][0][0],cube[1][0][1],relPos[2]),linterpolate(cube[1][10][0], cube[1][11][1],relPos[2]),relPos[1]),relPos[0]);
}

int main(int argc,char **args){
    SDL_Init(SDL_INIT_VIDEO);
    SDL_Surface *screen=SDL_SetVideoMode(512,512,32,SDL_SWSURFACE);
    srandom(SDL_GetTicks());  //If not on OSX/BSD, use srand()
    Uint32 *pixels;
    Uint32 grays[256];
    for(int x=0;x<256;x++){
        grays[x]=SDL_MapRGB(screen->format,x,x,x);
    }


    float grads[768];
    Uint8 perms[256];
    //First, generate the gradients and populate the permutation indexes.
    for(int x=0;x<256;x++){
        grads[3*x]=random();    //If not on OSX/BSD, use rand()
        grads[3*x+1]=random();
        grads[3*x+2]=random();
        normalize3(grads+3*x);

        perms[x]=x;
    }

    //Let's scramble those indexes!
    for(int x=0;x<256;x++){
        Uint8 temp=perms[x];
        Uint8 index=random();
        perms[x]=perms[index];
        perms[index]=temp;
    }

    printf("Permutation Indexes: ");
    for(int x=0;x<256;x++){
        printf("%hhu, ",perms[x]);
    }
    putchar('\n');

    Uint32 timer=SDL_GetTicks(),frameDelta;
    SDL_Event eventos;
    float zoom=-5.0;
    eventos.type=SDL_NOEVENT;
    while(eventos.type!=SDL_QUIT){
        SDL_PollEvent(&eventos);
        if(SDL_GetKeyState(NULL)[SDLK_UP]){
            zoom-=0.001*frameDelta;
        }
        else if(SDL_GetKeyState(NULL)[SDLK_DOWN]){
            zoom+=0.001*frameDelta;
        }
        float scale=expf(zoom);
        pixels=screen->pixels;
        float pos[3];
        pos[2]=SDL_GetTicks()/3000.0;
        for(int y=0;y<512;y++){
            pos[1]=y*scale;
            for(int x=0;x<512;x++){
                pos[0]=x*scale;
                float fracPos[3];
                fracPos[0]=pos[0];
                fracPos[1]=pos[1];
                fracPos[2]=pos[2];
                float color=noise3(fracPos,grads,perms);

                //Fractal sums of noise, if desired
                /*fracPos[0]*=2.0;
                fracPos[1]*=2.0;
                fracPos[2]*=2.0;
                color+=noise3(fracPos,grads,perms)/2.0;

                fracPos[0]*=2.0;
                fracPos[1]*=2.0;
                fracPos[2]*=2.0;
                color+=noise3(fracPos,grads,perms)/4.0;

                fracPos[0]*=2.0;
                fracPos[1]*=2.0;
                fracPos[2]*=2.0;
                color+=noise3(fracPos,grads,perms)/8.0;

                fracPos[0]*=2.0;
                fracPos[1]*=2.0;
                fracPos[2]*=2.0;
                color+=noise3(fracPos,grads,perms)/16.0;

                */

                *pixels++=grays[127+(Sint8)(256.0*color)];
            }
        }

        SDL_Flip(screen);
        frameDelta=SDL_GetTicks()-timer;
        printf("Running @ %.3f FPS!\n",1000.0/frameDelta);
        if(frameDelta<16){
            SDL_Delay(16-frameDelta);
        }
        timer=SDL_GetTicks();
    }

    return 0;
}

Usage: While running, press and hold Up or Down to magnify or shrink the noise grid.

Answer 1

I have finally found the problem: the gradient generator.

I was assuming that the random() function would pass its binary value to the grads[] array, covering the whole range of floating point numbers in the way. Unfortunately, that was not the case: its return value was being converted first into a float, then stored in the array. My biggest issue with this was that all generated vectors had positive member values.

This justified the block artifacts: there were lots of "hills" (high values) being generated next to each other, but no "valleys" (low values), and two adjacent hills would eventually clash and generate the lines along the integer values.

After realizing this, I tried to do some pointer juggling and storing the values directly in Uint32 form, but the values in the gradients became wacky (infs, NaNs, 1.0s and 0.0s all the way), so I came back to the original route and negated the numbers in the code itself.

This 7-liner solved the whole problem:

int y=random()&7;
if(y&1)
    grads[3*x]*=-1.0f;
if(y&2)
    grads[3*x+1]*=-1.0f;
if(y&4)
    grads[3*x+2]*=-1.0f;

Just place it before or after the normalizing function and it's done.

Now it looks like Perlin Noise:

And the fractal sum also looks a bit better:

@DiJuMx: I've seen the "Improving Noise" paper before, but hadn't realized how much the gradients would affect the noise appearance. Also, by trying to change the coordinate space from 0~256 to 0~1 resulted in the fractal sum not working anymore, and the resulting image had the same block artifacts.

Answer 2

It's a problem with the original implementation of Perlin noise.

He has a paper on it here

During the calculation of the gradient at integer coordinates, one or more of the vectors used will be 0, hence the overall gradient will be 0. As a result you get a grid of lines at the integer coordinates.

One way to fix this is to have the coordinate space go from 0 to 1 rather 0 to 512.

Another way is to implement the fix as describer in his paper.

Or finally, don't use the original Perlin Noise, instead swap to Simplex Noise which he also developed, paper here, and explanation here.