This code simulates hydraulic erosion on a greyscaled heightmap. Algorithm is taken from this paper (the optimized one): http://oddlabs.com/download/terrain_generation.pdf

Hydraulic erosion simulates changes to the terrain caused by flowing water dissolving material, transporting it and depositing it elsewhere.



Code:
//Hydraulic Erosion Filter
//recommended default values:
//iterations_ 		= 50
//threshold_		= 4096/heightmap size
var* TERRAIN_hydraulic_erosion_tmp=NULL;
///////////////////////////////
void TERRAIN_bmap_hydraulic_erosion_filter(BMAP* bmap_,var iterations_,var threshold_)
{	
	//get size
	double size_x_=bmap_width(bmap_);
	double size_y_=bmap_height(bmap_);
	
	//create data
	TERRAIN_hydraulic_erosion_tmp=(var*)sys_malloc(sizeof(var) * size_x_*size_y_);
	
	//fill data with pixel color
	int py; for(py=0;py<size_y_;py++)
	{
		int px; for(px=0;px<size_x_;px++)
		{
			COLOR color_;
			var format=bmap_lock(bmap_,0);var alpha_;
			var pixel=pixel_for_bmap(bmap_, px, py);
			pixel_to_vec(color_,alpha_,format,pixel);
			bmap_unlock(bmap_);
			
			TERRAIN_hydraulic_erosion_tmp[px*size_x_+py]=color_.red;
		}
	}
	
	int i; for(i=0;i<iterations_;i++)
	{
		int py; for(py=0;py<size_y_;py++)
		{
			int px; for(px=0;px<size_x_;px++)
			{
				var h_=TERRAIN_hydraulic_erosion_tmp[px*size_x_+py];
				
				var d[4]; COLOR color_tmp_; var h[4];
				
				h[0]=TERRAIN_hydraulic_erosion_tmp[clamp(px-1,0,size_x_-1)*size_x_+clamp(py,0,size_y_-1)];
				h[1]=TERRAIN_hydraulic_erosion_tmp[clamp(px,0,size_x_-1)*size_x_+clamp(py-1,0,size_y_-1)];
				h[2]=TERRAIN_hydraulic_erosion_tmp[clamp(px+1,0,size_x_-1)*size_x_+clamp(py,0,size_y_-1)];
				h[3]=TERRAIN_hydraulic_erosion_tmp[clamp(px,0,size_x_-1)*size_x_+clamp(py+1,0,size_y_-1)];

				var dtotal=0,dmax=0,l=0;
				int j; for(j=0;j<4;j++)
				{
					d[j]=h_-h[j];
					if(d[j]>dmax)
					{
						dmax=d[j];
						l=j;
					}
				}
				
				var hd;
				if(dmax>0 && dmax<=threshold_)
				{
					hd=0.5*dmax;
					h_=h_-hd;
					h[l]=h[l]+hd;
				}
				
				TERRAIN_hydraulic_erosion_tmp[px*size_x_+py]=h_;
				TERRAIN_hydraulic_erosion_tmp[clamp(px-1,0,size_x_-1)*size_x_+clamp(py,0,size_y_-1)]=h[0];
				TERRAIN_hydraulic_erosion_tmp[clamp(px,0,size_x_-1)*size_x_+clamp(py-1,0,size_y_-1)]=h[1];
				TERRAIN_hydraulic_erosion_tmp[clamp(px+1,0,size_x_-1)*size_x_+clamp(py,0,size_y_-1)]=h[2];
				TERRAIN_hydraulic_erosion_tmp[clamp(px,0,size_x_-1)*size_x_+clamp(py+1,0,size_y_-1)]=h[3];
			}
		}
	}
	
	//fill pixels with array data
	int py; for(py=0;py<size_y_;py++)
	{
		int px; for(px=0;px<size_x_;px++)
		{
			//set colors
			COLOR bmp_color_;
			bmp_color_.red=TERRAIN_hydraulic_erosion_tmp[px*size_x_+py];
			bmp_color_.green=TERRAIN_hydraulic_erosion_tmp[px*size_x_+py];
			bmp_color_.blue=TERRAIN_hydraulic_erosion_tmp[px*size_x_+py];
			
			//add to heightmap
			var format=bmap_lock(bmap_,0);
			var pixel=pixel_for_vec(bmp_color_,100,format);
			pixel_to_bmap(bmap_, px, py, pixel);
			bmap_unlock(bmap_);
		}
	}
	
	//free data
	sys_free(TERRAIN_hydraulic_erosion_tmp);
}