Problem with diamond-square algorithm

[kronikel]

I definitely agree it's the first few points that are messing it up.
And when you see a cone in the middle poking up, it's not because that point got set too high, it's because the next four points of the following square step got set too low.
I've considered a lot of ways to fix this, such as using a smoothing algorithm or a spline, or maybe a logarithm approach for how the roughness decays.
But I still can't get past how I see these "flag poles" even without randomness at all.
And what is the difference between my code and other ones?
When I look at other people's algorithms I basically see Height = Average of 4 surrounding points +- random
It feels like I'm missing something fundamental.

[knighty]

You are right. The spikes don't appear clearly in the 2D pictures but they are obvious in 3D. Mandelbrot's midpoint displacement algorithm doesn't have this issue. This quote from wikipedia suggests that the spikes issue is known:

The idea was first introduced by Fournier, Fussell and Carpenter at SIGGRAPH 1982. It was later analyzed by Gavin S. P. Miller in SIGGRAPH 1986 who described it as flawed.

@fractower - The randomness can't be the issue because when I take out the randomness completely, so it simply sets a height based on the average of the 4 corner heights, I get identical results.

That shouldn't happen because you are initializing the whole array height to the same value... Weird... Aren't you keeping the -r term when removing randomness? That will give a nice looking fractal surface anyway .

@David Makin: Thanks! I don't know how much slower FFT is but you can use arbitrary filters with it (more possibilities/variations).

[DarkBeam]

Already looked here? http://gameprogrammer.com/fractal.html#diamond

Luca

[kronikel]

I think it is considered flawed because the terrain produced isn't exactly 100% what is found in nature. But I've never been sure what he meant by that.

When I remove the randomness I also set the very middle point to be up, otherwise it would just create a completely flat surface.
But the middle point is the only one I touch and I had assumed I would get a smooth transition from the base up to the middle point, but I get some cones going down into the terrain.

@DarkBeam I used that exact method, or so I thought. I can't see what I'm doing different.

[David Makin]

@David Makin: Thanks! I don't know how much slower FFT is but you can use arbitrary filters with it (more possibilities/variations).

Yes, but if you look at the raw, basic Perlin method it allows an infinite number of ways of implementing it as the two components involved can be any type that fits the bill i.e. one noise function and one interpolation algorithm - these can essentially be any that you can conceive assuming they have the required number of dimensions (or can be adapted/expanded to the right number of dimensions).

See:

http://freespace.virgin.net/hugo.elias/models/m_perlin.htm

And remember that *any* "noise" will do and *any* "interpolation" will do, in fact really based on that as a simplified/generalised definition of Perlin Noise then it essentially includes the D-S algorithm.

[David Makin]

I think it is considered flawed because the terrain produced isn't exactly 100% what is found in nature. But I've never been sure what he meant by that.

When I remove the randomness I also set the very middle point to be up, otherwise it would just create a completely flat surface.
But the middle point is the only one I touch and I had assumed I would get a smooth transition from the base up to the middle point, but I get some cones going down into the terrain.

@DarkBeam I used that exact method, or so I thought. I can't see what I'm doing different.

I think I found the issue - here:

            int
                Size = Convert.ToInt32(System.Math.Pow(2, int.Parse(tbDetail.Text)) + 1),
                SideLength = Size - 1,
                d = 1025 / (Size - 1),
                HalfSide;
            Heights = new Point3D[Size, Size];


Size will be a power of two and SideLength a power of two minus 1, so if Size is 1024 then SideLength is 1023.
This will *not* work as SideLength itself *has* to be a power of 2, so you can either do this:

            int
                SideLength = Convert.ToInt32(System.Math.Pow(2, int.Parse(tbDetail.Text)) + 1),
                Size = SideLength + 1,
                d = 1025 / Size,
                HalfSide;
            Heights = new Point3D[Size, Size];

Or this:

            int
                SideLength = Convert.ToInt32(System.Math.Pow(2, int.Parse(tbDetail.Text)) + 1),
                Size = SideLength + 1,
                d = 1025 / Size,
                HalfSide;
            Heights = new Point3D[SideLength, SideLength];

but in this case also bitwise & the array-look-up index values with (SideLength-1).

[lycium]

[lycium]

rather evaluate it directly in a stateless, functional and parallel way

[kronikel]

I'm not sure what you mean.
With the way I have it right now SideLength is a power of 2 and Size is a power of 2 + 1.
So using "10" I have a Size of 1025 and a SideLength of 1024 at the start.

[cKleinhuis]

It isnt real because perlin or diamond squares wont peoduce caves or holes. . .

[David Makin]

It isnt real because perlin or diamond squares wont peoduce caves or holes. . .

Correct - but if you know a landscape generation algorithm that can return the height/heights per point in the horizontal plane *at runtime* as fast as using an appropriate version of D-S then please point me at it !!

(note I do mean arbitrary points not just an orderly set)

[kronikel]

I could be wrong, but I assumed things like that were due to weathering/ erosion/ etc.
I have seen people take terrain and apply another algorithm on top of it to make these though.

And about my problem, I applied a smoothing algorithm to my terrain, and while it did get rid of the cones, it also destroys the detail in the mountains and made everything look like smooth plastic.
So I changed it to only smooth out points that are higher/lower than the points around it by a certain amount, say 30 pixels.
Better results, but this is still just a patch and not nearly a fix for my problem.

[knighty]

It isnt real because perlin or diamond squares wont peoduce caves or holes. . .

A 3D version of diamond-square algorithm would generate caves . Though I haven't implemented it, that algorithm would require 3 stages:
-cube center point: average of 8 cube corners + random,
-faces center points: average of 4 corners and 2 adjacent cube centers + random.
-edges center points: average of 2 corners and 4 adjacent faces centers +random.
It could also be generelized to higher dimensions. For N dimension you'll need N stages.

And about my problem, I applied a smoothing algorithm to my terrain, and while it did get rid of the cones, it also destroys the detail in the mountains and made everything look like smooth plastic.
So I changed it to only smooth out points that are higher/lower than the points around it by a certain amount, say 30 pixels.
Better results, but this is still just a patch and not nearly a fix for my problem.

You could try using a filter that "detects the spikes". Here is a simple example:

Code:

filter(s){//s is the "strength" of the filter. A good value is around 10 if the heigtmap max height is around 1
   //edges are not processed here but they need to!
   for(l=1; l<N; l++)
      for(c=1; c<N; c++){
         G=tab[l][c]-tab[l-1][c]; G+=tab[l][c]-tab[l+1][c];
         G+=tab[l][c]-tab[l][c-1]; G+=tab[l][c]-tab[l][c+1];//"gradients" in the 4 directions. abs(G) will be higher on spike and "pits". 
         G=exp(-s*abs(G));//in order to have in in the range [-1, 1]
         tab[l][c]=G*tab[l][c]+0.25*(1-G)*(tab[l-1][c]+tab[l+1][c]+tab[l][c-1]+tab[l][c+1]);//simple filter
      }
}

I hope it doesn't smoothen the surface too much.

(In case someone is interested, I've attached some very simple evaldraw scripts about FFT based spectral synthesis noise (finally not that slow ); perlin noise (quite fast (EDIT: it is comparable to original perlin noise there was a bug with name collision) and GPU friendly (link) because it doesn't use any look up table), diamond-square and square-square algorithms)

[kronikel]

I made a little animation of the structure being created with no randomness at all, and the center point set up a little.
Hasn't really inspired a working idea for me yet, but maybe someone else can see what's going on.

[fractower]

That is a good animation Kronikel. It shows the wavelet basis functions for the square steps. There is a very similar basis for the diamond steps as well. The final surface is a linear combination of these basis functions. Your animation illustrates how the basis set is flawed. I have tried to address this by modifying the basis functions by performing a multigrid smoothing operation inside the main loop.  This required adding two new parameters. Erosion is a smoothing parameter and Erosion_delay is the number of main loops before smoothing is applied. The basis sets are for (Erosion_delay = 1, Erosion = .85) and (Erosion_delay = 2, Erosion = .55). I have also attached terrains made from the delay 2 basis and the original DS results.

The computation cost is ~1.5 times the original DS algorithm.

Code:

void model()  
{
  
  
  
  int
    Size = GRID,
    SideLength = Size - 1,
    SideLength2 = Size - 1,
    d = 1025 / (Size - 1),
    HalfSide, 
    HalfSide2,
    QSide2;
  //Heights = new Point3D[Size, Size];
  //float
  float r = r_init;
  float e = Erosion;
  //Roughness = float.Parse(RoughnessBox.Text);
  float hit;
  int grid_cnt = 0;

  srand(rand_init);

  for (int x = 0; x < Size; x++)
    for (int y = 0; y < Size; y++)
      Heights[x][y] = 0;

  
  while (SideLength2 >= 4)
    {
      HalfSide = SideLength / 2;
      HalfSide2 = SideLength2 / 2;
      QSide2 = HalfSide2 / 2;
      if(SideLength >=2){
	if(SideLength == 2) r = r / Roughness;
	for (int x = 0; x < Size - 1; x = x + SideLength)
	  {
	    for (int y = 0; y < Size - 1; y = y + SideLength)
	      {
		Heights[x + HalfSide][ y + HalfSide] =
		  (Heights[x][ y] +
		   Heights[x + SideLength][ y] +
		   Heights[x][ y + SideLength] +
		   Heights[x + SideLength][ y + SideLength]) / 4 + r - 2*r*(rand()/(1.0*RAND_MAX));
	      }
	  }
	
	for (int x = 0; x < Size - 1; x = x + SideLength)
	  {
	    for (int y = 0; y < Size - 1; y = y + SideLength)
	      {
		if (y != 0)
		  Heights[x + HalfSide][ y] = 
		    (Heights[x][ y] + Heights[x + SideLength][ y] + Heights[x + HalfSide][ y + HalfSide] + Heights[x + HalfSide][ y - HalfSide]) / 4 + r - 2*r*(rand()/(1.0*RAND_MAX)); 
		if (x != 0)
		  Heights[x][ y + HalfSide] = 
		    (Heights[x][ y] + Heights[x][ y + SideLength] + Heights[x + HalfSide][ y + HalfSide] + Heights[x - HalfSide][ y + HalfSide]) / 4 + r - 2*r*(rand()/(1.0*RAND_MAX));
	      }
	  }
      }
      
      if(grid_cnt >= Erosion_delay){
	for (int x = 0; x < Size - 1; x = x + SideLength2)
	  {
	    for (int y = 0; y < Size - 1; y = y + SideLength2)
	      {
		
		Heights[x + HalfSide2][ y + HalfSide2] = Heights[x + HalfSide2][ y + HalfSide2]*(1-e)+
		  (Heights[x + QSide2][ y + QSide2] +
		   Heights[x + SideLength2 - QSide2][ y + QSide2] +
		   Heights[x + QSide2][ y + SideLength2 - QSide2] +
		   Heights[x + SideLength2 - QSide2][ y + SideLength2 - QSide2]) * e / 4;
		
	      }
	  }
	
	for (int x = 0; x < Size - 1; x = x + SideLength2)
	  {
	    for (int y = 0; y < Size - 1; y = y + SideLength2)
	      {
		if (y != 0)
		  Heights[x + HalfSide2][ y] = Heights[x + HalfSide2][ y]*(1-e)+
		    (Heights[x + QSide2][ y] + Heights[x + SideLength2 - QSide2][ y] + 
		     Heights[x + HalfSide2][ y  + QSide2] + Heights[x + HalfSide2][ y - QSide2]) * e / 4; 
		if (x != 0)
		Heights[x][ y + HalfSide2] = Heights[x][ y + HalfSide2]*(1-e)+
		  (Heights[x][ y + QSide2] + Heights[x][ y + SideLength2 - QSide2] + 
		   Heights[x + QSide2][ y + HalfSide2] + Heights[x - QSide2][ y + HalfSide2]) * e / 4 ;
	      }
	  }
	SideLength2 = SideLength2 / 2;	
      }
      
      grid_cnt++;
      SideLength = SideLength / 2;
      r = r / Roughness;
    }
}