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 Yasha
(Posted 19 years ago)
!!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.
Open in online IDE
;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

Yasha
(Posted 18 years ago)
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