NewsForumsCodeLogsGallerySpecsSearch

Code archives/3D Graphics - Mesh/DSA heightmap

This code has been declared by its author to be Public Domain code.

Download source code

DSA heightmap by Yasha2008
!!Updated: See comments!!

Uses what I think is DSA. I hope it's useful to someone, it gives nice results for me. Does not attempt to include mesh saving or tiling features. 
;Another DSA heightmap alg.
;Seems to give nice results
;Tap space repeatedly to generate a tile
;LArrow and RArrow to turn camera

Graphics3D 800,600,32,2
SetBuffer BackBuffer()
WireFrame 1

centrecam=CreatePivot()
camera=CreateCamera(centrecam)


;Corner seed heights
y00=5
y0S=10
yS0=5
ySS=0

hconst#=5	;Height variance constant
size=32		;Grid size (should be a 2^n)
SeedRnd MilliSecs()


landscape=CreateMesh()
tile=CreateSurface(landscape)
EntityFX landscape,16		;so you can see it properly

iterations=Log(size)/Log(2)
iterated=0

Dim xpos(4)
Dim ypos#(4)
Dim zpos(4)
Dim dn(4)

PositionEntity centrecam,size/2,0,size/2
PositionEntity camera,16,15,-30
PointEntity camera,centrecam

Dim vert(size,size)
Dim poly(2*size,2*size)

cursor=CreateSphere()
EntityColor cursor,255,0,0
ScaleEntity cursor,0.5,0.5,0.5

For l=0 To size
	For r=0 To size
		vert(l,r)=AddVertex(tile,l,0,r)
	Next
Next

For l=1 To size
	For r=1 To size
		poly(l,r)=AddTriangle(tile,vert(l-1,r-1),vert(l-1,r),vert(l,r-1))
		poly(2*l,2*r)=AddTriangle(tile,vert(l-1,r),vert(l,r),vert(l,r-1))
	Next
Next

VertexCoords tile,vert(0,0),VertexX(tile,vert(0,0)),y00,VertexZ(tile,vert(0,0))
VertexCoords tile,vert(0,size),VertexX(tile,vert(0,size)),y0S,VertexZ(tile,vert(0,size))
VertexCoords tile,vert(size,0),VertexX(tile,vert(size,0)),yS0,VertexZ(tile,vert(size,0))
VertexCoords tile,vert(size,size),VertexX(tile,vert(size,size)),ySS,VertexZ(tile,vert(size,size))

While Not KeyDown(1)
	
	TurnEntity centrecam,0,KeyDown(205)-KeyDown(203),0
	;If KeyHit(28)=1 Then wired=1-wired
	;WireFrame wired
	
	If KeyHit(57)=1 And iterated<iterations
		For n=1 To 2^iterated
			For itn=1 To 2^iterated
				xpos(0)=(size/(2^(iterated+1)))+(n-1)*(size/(2^iterated))					;get location...
				zpos(0)=(size/(2^(iterated+1)))+(itn-1)*(size/(2^iterated))					;x1z1 etc give the square locations
				xpos(1)=xpos(0):xpos(3)=xpos(0)
				zpos(2)=zpos(0):zpos(4)=zpos(0)												;x4z1, x2z1, x2z3, x4z3 give the previous
				xpos(2)=n*(size/(2^iterated)):zpos(1)=itn*(size/(2^iterated))				;diamond of which x0z0 is centre
				xpos(4)=(n-1)*(size/(2^iterated)):zpos(3)=(itn-1)*(size/(2^iterated))
				
				ypos(0)=(VertexY(tile,vert(xpos(4),zpos(1)))+VertexY(tile,vert(xpos(2),zpos(1)))+VertexY(tile,vert(xpos(2),zpos(3)))+VertexY(tile,vert(xpos(4),zpos(3))))/4
				ypos(0)=ypos(0)+hvar(hconst,iterated)
				VertexCoords(tile,vert(xpos(0),zpos(0)),xpos(0),ypos(0),zpos(0))			;centres
				
				For q=1 To 4
					dn(q)=0
					ypos(q)=0
				Next
				
				
				If zpos(1)+(zpos(1)-zpos(0))<=size Then ypos(1)=ypos(1)+VertexY(tile,vert(xpos(1),zpos(1)+(zpos(1)-zpos(0)))):dn(1)=dn(1)+1
				ypos(1)=ypos(1) + VertexY(tile,vert(xpos(4),zpos(1))) + VertexY(tile,vert(xpos(2),zpos(1))) + ypos(0)  :dn(1)=dn(1)+3
				ypos(1)=ypos(1)/dn(1);:dn(1)=0
				ypos(1)=ypos(1)+hvar(hconst,iterated)
				
				If xpos(2)+(xpos(2)-xpos(0))<=size Then ypos(2)=ypos(2)+VertexY(tile,vert(xpos(2)+(xpos(2)-xpos(0)),zpos(2))):dn(2)=dn(2)+1
				ypos(2)=ypos(2)+VertexY(tile,vert(xpos(2),zpos(1)))+VertexY(tile,vert(xpos(2),zpos(3)))+ypos(0):dn(2)=dn(2)+3
				ypos(2)=ypos(2)/dn(2);:dn(2)=0
				ypos(2)=ypos(2)+hvar(hconst,iterated)
				
				If zpos(3)-(zpos(3)-zpos(0))>=0 Then ypos(3)=ypos(3)+VertexY(tile,vert(xpos(3),zpos(3)-(zpos(3)-zpos(0)))):dn(3)=dn(3)+1
				ypos(3)=ypos(3)+VertexY(tile,vert(xpos(4),zpos(3)))+VertexY(tile,vert(xpos(2),zpos(3)))+ypos(0):dn(3)=dn(3)+3
				ypos(3)=ypos(3)/dn(3);:dn(3)=0
				ypos(3)=ypos(3)+hvar(hconst,iterated)
				
				If xpos(4)-(xpos(4)-xpos(0))>=0 Then ypos(4)=ypos(4)+VertexY(tile,vert(xpos(4)-(xpos(4)-xpos(0)),zpos(4))):dn(4)=dn(4)+1
				ypos(4)=ypos(4)+VertexY(tile,vert(xpos(4),zpos(1)))+VertexY(tile,vert(xpos(4),zpos(3)))+ypos(0):dn(4)=dn(4)+3
				ypos(4)=ypos(4)/dn(4);:dn(4)=0
				ypos(4)=ypos(4)+hvar(hconst,iterated)
				
				VertexCoords(tile,vert(xpos(1),zpos(1)),xpos(1),ypos(1),zpos(1))			;squares
				VertexCoords(tile,vert(xpos(2),zpos(2)),xpos(2),ypos(2),zpos(2))
				VertexCoords(tile,vert(xpos(3),zpos(3)),xpos(3),ypos(3),zpos(3))
				VertexCoords(tile,vert(xpos(4),zpos(4)),xpos(4),ypos(4),zpos(4))
			Next
		Next
		iterated=iterated+1
	EndIf
	
	RenderWorld

	Text 0,0,iterated
	Text 0,20,2^iterated
	Text 0,40,n
	For q=0 To 4
		Text 50,20*q,ypos(q)
		Text 200,20*q,dn(q)
	Next

Flip
Wend
End

Function hvar#(k#,itn)

	range#=(k/2^itn)
	variance#=Rnd(range*-1,range)
	Return variance

End Function

Comments

Yasha2009
The above version is clearer if you want to understand the algorithm, but not very portable. Here's the same thing inside a function, no globals or Dims required! Returns an image rather than a mesh, so you can use it in 2D too.

Function DSA(fx,fz,hconst#,subk=0,jagk#=0.5,density#=1,maxh#=1,seed=-1)		;Final width, final height, height variance, subfeatures, softness
																		;density, alpha and RNDgen seed (so you can repeat a specific pattern)
Local iterations,iterated#,jagv#,hipoint#,lopoint#,dn,size,max,sig,rs,scale#,t#
Local xpos[4],ypos[4],zpos[4],dpos[4],dput[1],dd[1]
Local mx=2^Ceil(Log(fx)/Log(2))
Local mz=2^Ceil(Log(fz)/Log(2))

If mx>mz Then size=mx:Else size=mz
iterations=Log(size)/Log(2)
max=size+1:map=CreateBank(max*max*4)

If seed>-1 Then rs=RndSeed():SeedRnd seed	;So that other functions stay predictable

For iterated=0 To iterations-1
	jagv=jagk*((iterated/iterations)^2)
	For n=1 To 2^iterated
		For itn=1 To 2^iterated
			xpos[0]=(size/(2^(iterated+1)))+(n-1)*(size/(2^iterated))					;get location...
			zpos[0]=(size/(2^(iterated+1)))+(itn-1)*(size/(2^iterated))					;x1z1 etc give the square locations
			xpos[1]=xpos[0]:xpos[3]=xpos[0]
			zpos[2]=zpos[0]:zpos[4]=zpos[0]												;x4z1, x2z1, x2z3, x4z3 give the previous
			xpos[2]=n*(size/(2^iterated)):zpos[1]=itn*(size/(2^iterated))				;diamond of which x0z0 is centre
			xpos[4]=(n-1)*(size/(2^iterated)):zpos[3]=(itn-1)*(size/(2^iterated))
			
			dput[0]=0:dpos[2]=max*xpos[4]+zpos[3]:dpos[3]=max*xpos[2]+zpos[1]:dd[0]=2			;Extended diagonals for quadratic/cubic smoothing
			If xpos[4]-((xpos[0]-xpos[4])*2)>=0 And zpos[3]-((zpos[0]-zpos[3])*2)>=0 Then dpos[1]=max*(xpos[4]-((xpos[0]-xpos[4])*2))+(zpos[3]-((zpos[0]-zpos[3])*2)):dd[0]=3
			If xpos[2]+((xpos[2]-xpos[0])*2)<=size And zpos[1]+((zpos[1]-zpos[0])*2)<=size Then dpos[4]=max*(xpos[2]+((xpos[2]-xpos[0])*2))+(zpos[1]+((zpos[1]-zpos[0])*2)):dd[0]=dd[0]+2
			If dd[0]>2
				If dd[0]=3 Then dput[0]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[0]=4 Then dput[0]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[0]=5 Then dput[0]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			
			dpos[2]=max*xpos[2]+zpos[3]:dpos[3]=max*xpos[4]+zpos[1]:dd[1]=2			;Opposite extended diagonals for quadratic/cubic smoothing
			If xpos[2]+((xpos[2]-xpos[0])*2)<=size And zpos[3]-((zpos[0]-zpos[3])*2)>=0 Then dpos[1]=max*(xpos[2]+((xpos[2]-xpos[0])*2))+(zpos[3]-((zpos[0]-zpos[3])*2)):dd[1]=3
			If xpos[4]-((xpos[0]-xpos[4])*2)>=0 And zpos[1]+((zpos[1]-zpos[0])*2)<=size Then dpos[4]=max*(xpos[4]-((xpos[0]-xpos[4])*2))+(zpos[1]+((zpos[1]-zpos[0])*2)):dd[1]=dd[1]+2
			If dd[1]>2
				If dd[1]=3 Then dput[1]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[1]=4 Then dput[1]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[1]=5 Then dput[1]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			
			bleh1=xpos[2]
			bleh2=zpos[3]
			
			ypos[0]=(PeekHeight(map,max*xpos[4]+zpos[1])+PeekHeight(map,max*xpos[2]+zpos[1])+PeekHeight(map,max*xpos[2]+zpos[3])+PeekHeight(map,max*xpos[4]+zpos[3]))/4
			ypos[0]=ypos[0]+hvar(hconst,iterated,subk)
			If dd[0]+dd[1]>4
				If dd[0]>2 And dd[1]>2
					ypos[0]=linpol(ypos[0],(dput[0]+dput[1])/2,jagv)
				Else
					dput[0]=linpol(ypos[0],dput[0]+dput[1],jagv)
					ypos[0]=linpol(ypos[0],dput[0],jagv)
				EndIf
			EndIf
			PokeHeight map,max*xpos[0]+zpos[0],ypos[0]
			If ypos[0]<lopoint Then lopoint=ypos[0]
			If ypos[0]>hipoint Then hipoint=ypos[0]
					
			dn=0:dd[0]=0:dd[1]=0
			dput[0]=0:dput[1]=0
			ypos[1]=0:dpos[1]=0
			ypos[2]=0:dpos[2]=0
			ypos[3]=0:dpos[3]=0
			ypos[4]=0:dpos[4]=0
			
			
			If zpos[1]+(zpos[1]-zpos[0])<=size
				ypos[1]=ypos[1]+PeekHeight(map,max*xpos[1]+(zpos[1]+(zpos[1]-zpos[0])))
				dn=1
				dpos[3]=max*xpos[1]+(zpos[1]+(zpos[1]-zpos[0]))
				dd[0]=1
			EndIf
			dpos[2]=max*xpos[1]+zpos[0]:dd[0]=dd[0]+1
			ypos[1]=ypos[1]+PeekHeight(map,max*xpos[4]+zpos[1])+PeekHeight(map,max*xpos[2]+zpos[1])+ypos[0]:dn=dn+3
			ypos[1]=ypos[1]/dn:dn=0
			ypos[1]=ypos[1]+hvar(hconst,iterated,subk)
			If dd[0]>1
				If zpos[0]-(2*(zpos[1]-zpos[0]))>=0 Then dpos[1]=max*xpos[1]+(zpos[0]-(2*(zpos[1]-zpos[0]))):dd[0]=dd[0]+1
				If zpos[1]+(3*(zpos[1]-zpos[0]))<=size Then dpos[4]=max*xpos[1]+(zpos[1]+(3*(zpos[1]-zpos[0]))):dd[0]=dd[0]+2
			EndIf
			If dd[0]>2
				If dd[0]=3 Then dput[0]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[0]=4 Then dput[0]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[0]=5 Then dput[0]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			dpos[2]=max*xpos[4]+zpos[1]:dpos[3]=max*xpos[2]+zpos[1]:dd[1]=2
			If xpos[4]-(2*(xpos[1]-xpos[4]))>=0 Then dpos[1]=max*(xpos[4]-(2*(xpos[1]-xpos[4])))+zpos[1]:dd[1]=dd[1]+1
			If xpos[2]+(2*(xpos[2]-xpos[1]))<=size Then dpos[4]=max*(xpos[2]+(2*(xpos[2]-xpos[1])))+zpos[1]::dd[1]=dd[1]+2
			If dd[1]>2
				If dd[1]=3 Then dput[1]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[1]=4 Then dput[1]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[1]=5 Then dput[1]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			If dd[0]+dd[1]>3
				If dd[0]>2 And dd[1]>2
					ypos[1]=linpol(ypos[1],(dput[0]+dput[1])/2,jagv)
				Else
					dput[0]=linpol(ypos[1],dput[0]+dput[1],jagv)
					ypos[1]=linpol(ypos[1],dput[0],jagv)
				EndIf
			EndIf
			dd[0]=0:dd[1]=0
			dput[0]=0:dput[1]=0
			If ypos[1]<lopoint Then lopoint=ypos[1]
			If ypos[1]>hipoint Then hipoint=ypos[1]
			
			If xpos[2]+(xpos[2]-xpos[0])<=size
				ypos[2]=ypos[2]+PeekHeight(map,max*(xpos[2]+(xpos[2]-xpos[0]))+zpos[2])
				dn=1
				dpos[3]=max*(xpos[2]+(xpos[2]-xpos[0]))+zpos[2]
				dd[0]=1
			EndIf
			dpos[2]=max*xpos[0]+zpos[2]:dd[0]=dd[0]+1
			ypos[2]=ypos[2]+PeekHeight(map,max*xpos[2]+zpos[1])+PeekHeight(map,max*xpos[2]+zpos[3])+ypos[0]:dn=dn+3
			ypos[2]=ypos[2]/dn:dn=0
			ypos[2]=ypos[2]+hvar(hconst,iterated,subk)
			If dd[0]>1
				If xpos[0]-(2*(xpos[2]-xpos[0]))>=0 Then dpos[1]=max*(xpos[0]-(2*(xpos[2]-xpos[0])))+zpos[2]:dd[0]=dd[0]+1
				If xpos[2]+(3*(xpos[2]-xpos[0]))<=size Then dpos[4]=max*(xpos[2]+(3*(xpos[2]-xpos[0])))+zpos[2]:dd[0]=dd[0]+2
			EndIf
			If dd[0]>2
				If dd[0]=3 Then dput[0]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[0]=4 Then dput[0]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[0]=5 Then dput[0]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			dpos[2]=max*xpos[2]+zpos[3]:dpos[3]=max*xpos[2]+zpos[1]:dd[1]=2
			If zpos[3]-(2*(zpos[2]-zpos[3]))>=0 Then dpos[1]=max*xpos[2]+(zpos[3]-(2*(zpos[2]-zpos[3]))):dd[1]=dd[1]+1
			If zpos[1]+(2*(zpos[1]-zpos[2]))<=size Then dpos[4]=max*xpos[2]+(zpos[1]+(2*(zpos[1]-zpos[2]))):dd[1]=dd[1]+2
			If dd[1]>2
				If dd[1]=3 Then dput[1]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[1]=4 Then dput[1]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[1]=5 Then dput[1]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			If dd[0]+dd[1]>3
				If dd[0]>2 And dd[1]>2
					ypos[2]=linpol(ypos[2],(dput[0]+dput[1])/2,jagv)
				Else
					dput[0]=linpol(ypos[2],dput[0]+dput[1],jagv)
					ypos[2]=linpol(ypos[2],dput[0],jagv)
				EndIf
			EndIf
			dd[0]=0:dd[1]=0
			dput[0]=0:dput[1]=0
			If ypos[2]<lopoint Then lopoint=ypos[2]
			If ypos[2]>hipoint Then hipoint=ypos[2]
			
			If zpos[3]-(zpos[0]-zpos[3])>=0
				ypos[3]=ypos[3]+PeekHeight(map,max*xpos[3]+(zpos[3]-(zpos[0]-zpos[3])))
				dn=1
				dpos[2]=max*xpos[3]+(zpos[3]-(zpos[0]-zpos[3]))
				dd[0]=1
			EndIf
			dpos[3]=max*xpos[3]+zpos[0]:dd[0]=dd[0]+1
			ypos[3]=ypos[3]+PeekHeight(map,max*xpos[4]+zpos[3])+PeekHeight(map,max*xpos[2]+zpos[3])+ypos[0]:dn=dn+3
			ypos[3]=ypos[3]/dn:dn=0
			ypos[3]=ypos[3]+hvar(hconst,iterated,subk)
			If dd[0]>1
				If zpos[3]-(3*(zpos[0]-zpos[3]))>=0 Then dpos[1]=max*xpos[3]+(zpos[3]-(3*(zpos[0]-zpos[3]))):dd[0]=dd[0]+1
				If zpos[0]+(2*(zpos[0]-zpos[3]))<=size Then dpos[4]=max*xpos[3]+(zpos[0]+(2*(zpos[0]-zpos[3]))):dd[0]=dd[0]+2
			EndIf
			If dd[0]>2
				If dd[0]=3 Then dput[0]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[0]=4 Then dput[0]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[0]=5 Then dput[0]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			dpos[2]=max*xpos[4]+zpos[3]:dpos[3]=max*xpos[2]+zpos[3]:dd[1]=2
			If xpos[4]-(xpos[2]-xpos[4])>=0 Then dpos[1]=max*(xpos[4]-(xpos[2]-xpos[4]))+zpos[3]:dd[1]=dd[1]+1
			If xpos[2]+(xpos[2]-xpos[4])<=size Then dpos[4]=max*(xpos[2]+(xpos[2]-xpos[4]))+zpos[3]:dd[1]=dd[1]+2
			If dd[1]>2
				If dd[1]=3 Then dput[1]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[1]=4 Then dput[1]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[1]=5 Then dput[1]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			If dd[0]+dd[1]>3
				If dd[0]>2 And dd[1]>2
					ypos[3]=linpol(ypos[3],(dput[0]+dput[1])/2,jagv)
				Else
					dput[0]=linpol(ypos[3],dput[0]+dput[1],jagv)
					ypos[3]=linpol(ypos[3],dput[0],jagv)
				EndIf
			EndIf
			dd[0]=0:dd[1]=0
			dput[0]=0:dput[1]=0
			If ypos[3]<lopoint Then lopoint=ypos[3]
			If ypos[3]>hipoint Then hipoint=ypos[3]
			
			If xpos[4]-(xpos[0]-xpos[4])>=0
				ypos[4]=ypos[4]+PeekHeight(map,max*(xpos[4]-(xpos[0]-xpos[4]))+zpos[4])
				dn=1
				dpos[2]=max*(xpos[4]-(xpos[0]-xpos[4]))+zpos[4]
				dd[0]=1
			EndIf
			dpos[3]=max*xpos[0]+zpos[4]:dd[0]=dd[0]+1
			ypos[4]=ypos[4]+PeekHeight(map,max*xpos[4]+zpos[1])+PeekHeight(map,max*xpos[4]+zpos[3])+ypos[0]:dn=dn+3
			ypos[4]=ypos[4]/dn:dn=0
			ypos[4]=ypos[4]+hvar(hconst,iterated,subk)
			If dd[0]>1
				If xpos[4]-(3*(xpos[0]-xpos[4]))>=0 Then dpos[1]=max*(xpos[4]-(3*(xpos[0]-xpos[4])))+zpos[4]:dd[0]=dd[0]+1
				If xpos[0]+(xpos[2]-xpos[4])<=size Then dpos[4]=max*(xpos[0]+(xpos[2]-xpos[4]))+zpos[4]:dd[0]=dd[0]+2
			EndIf
			If dd[0]>2
				If dd[0]=3 Then dput[0]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[0]=4 Then dput[0]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[0]=5 Then dput[0]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			dpos[2]=max*xpos[4]+zpos[3]:dpos[3]=max*xpos[4]+zpos[1]:dd[1]=2
			If zpos[3]-(zpos[1]-zpos[3])>=0 Then dpos[1]=max*xpos[4]+(zpos[3]-(zpos[1]-zpos[3])):dd[1]=dd[1]+1
			If zpos[1]+(zpos[1]-zpos[3])<=size Then dpos[4]=max*xpos[4]+(zpos[1]+(zpos[1]-zpos[3])):dd[1]=dd[1]+2
			If dd[1]>2
				If dd[1]=3 Then dput[1]=quadpol(dpos[1],dpos[2],dpos[3],max*xpos[0]+zpos[0],map,max)
				If dd[1]=4 Then dput[1]=quadpol(dpos[2],dpos[3],dpos[4],max*xpos[0]+zpos[0],map,max)
				If dd[1]=5 Then dput[1]=cubepol(dpos[1],dpos[2],dpos[3],dpos[4],map)
			EndIf
			If dd[0]+dd[1]>3
				If dd[0]>2 And dd[1]>2
					ypos[4]=linpol(ypos[4],(dput[0]+dput[1])/2,jagv)
				Else
					dput[0]=linpol(ypos[4],dput[0]+dput[1],jagv)
					ypos[4]=linpol(ypos[4],dput[0],jagv)
				EndIf
			EndIf
			dd[0]=0:dd[1]=0
			dput[0]=0:dput[1]=0
			If ypos[4]<lopoint Then lopoint=ypos[4]
			If ypos[4]>hipoint Then hipoint=ypos[4]
			
			PokeHeight(map,max*xpos[1]+zpos[1],ypos[1])
			PokeHeight(map,max*xpos[2]+zpos[2],ypos[2])
			PokeHeight(map,max*xpos[3]+zpos[3],ypos[3])
			PokeHeight(map,max*xpos[4]+zpos[4],ypos[4])
		Next
	Next
Next

scale=255.0/(hipoint-lopoint)

hipoint=hipoint-lopoint
For x=0 To size
For y=0 To size
	t=PeekHeight(map,max*x+y)
	t=t-lopoint
	t=t-(1-density)*hipoint
	If t<0 Then t=0
	PokeHeight(map,max*x+y,t)
Next
Next

hipoint=hipoint*density
scale=(255.0/hipoint)*maxh

dsamap=CreateImage(fx,fz)
LockBuffer ImageBuffer(dsamap)

For x=0 To fx-1
	For y=0 To fz-1
		hfinal#=PeekHeight(map,(max*x+y))*scale
		WritePixelFast(x,y,$ff000000 Or (hfinal Shl 16) Or (hfinal Shl 8) Or hfinal,ImageBuffer(dsamap))
	Next
Next

UnlockBuffer ImageBuffer(dsamap)
FreeBank map

If seed>-1 Then SeedRnd rs		;Restore to where it was before...
Return dsamap

End Function

Function PeekHeight(bank,offset)
	Return PeekFloat(bank,offset*4)
End Function

Function PokeHeight(bank,offset,value#)
	PokeFloat bank,offset*4,value#
End Function

Function hvar#(k#,itn,f=0)
	Local range#
	itn=itn-f
	If itn<1 Then itn=1
	range=k/(2^itn)
	Return Rnd(range*-1,range)
End Function

Function linpol#(a#,b#,x#)				;Linear interpolation
	Return (a*(1-x))+(b*x)
End Function

Function quadpol#(d,e,f,xp,bank,max)	;Quadratic interpolation
	Local a#,b#,c#,x#
		
	If e Mod max=xp Mod max
		x=0.5*Sgn(xp/max-e/max)
	Else
		x=0.5*Sgn((xp Mod max)-(e Mod max))
	EndIf
	
	c#=PeekHeight(bank,e)
	b#=(PeekHeight(bank,f)-PeekHeight(bank,d))/2
	a#=(PeekHeight(bank,d)-c)+b
	Return (a*(x*x))+(b*x)+c
End Function

Function cubepol#(v0,v1,v2,v3,bank)		;Cubic interpolation
	Local P#=(PeekHeight(bank,v3)-PeekHeight(bank,v2))-(PeekHeight(bank,v0)-PeekHeight(bank,v1))
	Local Q#=(PeekHeight(bank,v0)-PeekHeight(bank,v1))-P
	Local R#=PeekHeight(bank,v2)-PeekHeight(bank,v0)
	Local S#=PeekHeight(bank,v1)
	Return (P*0.125)+(Q*0.25)+(R*0.5)+S
End Function



Code Archives Forum