[4.60] Telluric Planet Generator

Share your advanced PureBasic knowledge/code with the community.
Kelebrindae
Enthusiast
Enthusiast
Posts: 151
Joined: Tue Apr 01, 2008 3:23 pm

[4.60] Telluric Planet Generator

Post by Kelebrindae »

Hi,

Here is a updated version of my "3D Planet generator". It's far from perfect and the code is quite messy, but the funny thing is that it allows to change the mesh LOD by tesselating/de-tesselating the mesh at runtime (by pressing R / T).
Image
(screenshot by applePi :wink: )
Have fun!

Code: Select all

; Author: Kelebrindae
; Date: december 27 2011
; PB version: v4.61

; -------------------------------------------------------------------------------------
; Purpose: Generate a random telluric planet
; -------------------------------------------------------------------------------------
; - To make make mapping easier and avoid distorsion to the poles, the mesh is actually
;   made of 6 planes that form a cube => all we need is 6 square heightmaps for the
;   relief, and 6 square textures.
; -------------------------------------------------------------------------------------
; Known bugs: 
; - Visible seams between heightmaps
; - Visible seams between submeshes
; -------------------------------------------------------------------------------------

;- ---------- Controls ----------
;- Arrow keys to rotate
;- W : wireframe display
;- R/T : decrease/increase tesselation
;- Mouse wheel: zoom
;- ------------------------------

;- Constants
#PLANETSEED=123456789 ; change this to create a different planet
#WATERLEVEL=127
#RELIEFEXAGERATION=0.2 ; Change this to increase/decrease relief height
#USEAPIDRAWING = #False
#GRIDSIZE=256 ; heightmaps and textures size

; Sides of the cube
Enumeration
  #TOPSIDE
  #BOTTOMSIDE
  #FRONTSIDE
  #BACKSIDE
  #LEFTSIDE
  #RIGHTSIDE
EndEnumeration

; Side of an image (to blend seams)
#HEIGHTMAP_TOP = 1
#HEIGHTMAP_BOTTOM = 2
#HEIGHTMAP_LEFT = 4
#HEIGHTMAP_RIGHT = 8


; Texture type
Enumeration
  #MAPTYPE_TERRAIN
  #MAPTYPE_ALTITUDE
  #MAPTYPE_SLOPE
  #MAPTYPE_TEMPERATURE
EndEnumeration


;- Data structures
Structure Vector3_struct ; used in vector maths (normalization)
  x.f
  y.f
  z.f
EndStructure

Structure VoronoiPoint2D_struct ; used to create the Voronoi Diagram (tectonic plates)
  x.i
  y.i
  dist.f
EndStructure

Structure geostat_struct ; used to store infos about terrains
  latitude.f
  slope.i
  temperature.f
  height.i
  probaterrain.f[10]
  isWater.b
  colour.i
EndStructure

Structure terrain_struct ; used to store specifications for each terrain type
  name.s
  heightAv.f
  heightVt.f
  slopeAv.f
  slopeVt.f
  temperatureAv.f
  temperatureVt.f
  colour.i
EndStructure

Structure vertex_struct
  x.f
  y.f
  z.f
  nx.f
  ny.f
  nz.f
  colour.i
  u.f
  v.f
EndStructure

Structure triangle_struct
  v1.i
  v2.i
  v3.i
 
  *ptrV1.vertex_struct
  *ptrV2.vertex_struct
  *ptrV3.vertex_struct
EndStructure

; API drawing
CompilerIf #USEAPIDRAWING = #True
  Structure myBITMAPINFO
    bmiHeader.BITMAPINFOHEADER
    bmiColors.RGBQUAD[1]
  EndStructure
CompilerEndIf 


;- Global definitions
Global CameraMode.b
Global anglex.f = 0,angley.f = 0


Global Dim diamond.b(#GRIDSIZE,#GRIDSIZE)
Global Dim voronoi.b(#GRIDSIZE,#GRIDSIZE)

; There are 6 heightmaps and geostats arrays: one for each side of the cube
Global Dim heightmap.i(6,#GRIDSIZE,#GRIDSIZE)
Global Dim newhmap.i(6,#GRIDSIZE,#GRIDSIZE)
Global Dim geostat.geostat_struct(6,#GRIDSIZE,#GRIDSIZE)

;- Terrain types
Global Dim terrain.terrain_struct(10)

; Temperatures for poles and equator (Celsius)
Global tpol.f = -15,teq.f = 50,name.s

; Vertices which are on the heightmaps' borders
Structure borderVert_struct
  *ptrVert.vertex_struct
  u.f
  v.f
EndStructure

; Mesh info for each side of the cube
Structure planetSide_struct
  numTexture.i
  numMaterial.i
  
  List topBorder.borderVert_struct()
  List bottomBorder.borderVert_struct()
  List leftBorder.borderVert_struct()
  List rightBorder.borderVert_struct()
  
  List vertex.vertex_struct()
  List triangle.triangle_struct()
EndStructure
Global Dim planetSide.planetSide_struct(6)

; Planet entity
Global planetEntity.i, tesselation.i = 5


EnableExplicit

;************************************************************************************
;-                                 ---- Macros ----
;************************************************************************************

Global SQR_sqrhalf.f,*SQR_sqrptr,SQR_sqrInt.i
; Fast inverse square root, Romero style
Macro FAST_INVSQR(x,result)
  result = x
  SQR_sqrhalf = 0.5 * result
  *SQR_sqrptr = @result
  SQR_sqrInt  = PeekI(*SQR_sqrptr)
  
  ; The magic number is different in 64 bits
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x64
    PokeI(*SQR_sqrptr,$5FE6EB50C7AA19F9 - ( SQR_sqrInt >> 1))
  CompilerElse ; all the others are 32 bits ?
    PokeI(*SQR_sqrptr,$5F375A86 - ( SQR_sqrInt >> 1)) ; or $5f3759df
  CompilerEndIf
  
  result=result*(1.5 - SQR_sqrhalf * result * result)
  result=result*(1.5 - SQR_sqrhalf * result * result) ; twice for greater precision
  result=result*(1.5 - SQR_sqrhalf * result * result) ; thrice is even better (but you can ditch it to squeeze a few extra cycles from your code)
EndMacro

; Vector3 normalization
Global VEC3_length.f
Macro VEC3_NORMALIZE(vector)
  FAST_INVSQR((vector\x*vector\x) + (vector\y*vector\y)+(vector\z*vector\z),VEC3_length)
  vector\x * VEC3_length
  vector\y * VEC3_length
  vector\z * VEC3_length
EndMacro 

Macro VEC3_NORMALIZETOLENGTH(vector,magnitude)
  FAST_INVSQR((vector\x*vector\x) + (vector\y*vector\y)+(vector\z*vector\z),VEC3_length)
  VEC3_length*(magnitude)
  vector\x * VEC3_length
  vector\y * VEC3_length
  vector\z * VEC3_length
EndMacro 

; Commpute a 3D point between two others (for tesselation)
Macro MIDDLEPOINT(result,pt1,pt2)
  result\x = pt1\x + (pt2\x - pt1\x) / 2
  result\y = pt1\y + (pt2\y - pt1\y) / 2
  result\z = pt1\z + (pt2\z - pt1\z) / 2
 
  result\u = pt1\u + (pt2\u - pt1\u) / 2
  result\v = pt1\v + (pt2\v - pt1\v) / 2
EndMacro

; Add a point in the list, only if it doesn't already exists. If it exists, return the existing point
Macro CREATENEWVERTEX(refPt,numNewPt)
  
  ; This uses the concatenated X,Y,Z coords as a map key to check for duplicate vertices (yeah, I know, it's ugly)
  ; ("Strf(...,2)" assumes that vertex coords within the 0.01 threshold are identical)
  numNewPt = antidup(StrF(refPt\x,2)+StrF(refPt\y,2)+StrF(refPt\z,2))-1
  If numNewPt = -1
    AddElement(newVertex())
    newVertex()\x = refPt\x
    newVertex()\y = refPt\y
    newVertex()\z = refPt\z
    newVertex()\nx = refPt\nx
    newVertex()\ny = refPt\ny
    newVertex()\nz = refPt\nz
    newVertex()\u = refPt\u
    newVertex()\v = refPt\v
    newVertex()\colour = refPt\colour
   
    numNewPt = ListIndex(newVertex())
    AddMapElement(antidup(),StrF(refPt\x,2)+StrF(refPt\y,2)+StrF(refPt\z,2),#PB_Map_NoElementCheck)
    antidup() = numNewPt+1
  EndIf
    
EndMacro

; Weld the vertices between two heightmaps.
; Ex: weld the bottom of the "Top" heightmap to the top of the "Front" heightmap
Macro WELDBORDERVERTICES(side1,border1,coord1,side1Sort,side2,border2,coord2)
  SortStructuredList(planetSide(side1)\border1#Border(),side1Sort,OffsetOf(borderVert_struct\coord1),#PB_Sort_Float)
  SortStructuredList(planetSide(side2)\border2#Border(),#PB_Sort_Ascending,OffsetOf(borderVert_struct\coord2),#PB_Sort_Float)
  
  FirstElement(planetSide(side2)\border2#Border())
  ForEach planetSide(side1)\border1#Border()    
    planetSide(side2)\border2#Border()\ptrVert\x = planetSide(side1)\border1#Border()\ptrVert\x
    planetSide(side2)\border2#Border()\ptrVert\y = planetSide(side1)\border1#Border()\ptrVert\y
    planetSide(side2)\border2#Border()\ptrVert\z = planetSide(side1)\border1#Border()\ptrVert\z

    NextElement(planetSide(side2)\border2#Border())
  Next planetSide(side1)\border1#Border()

EndMacro

;************************************************************************************
;-                                 ---- Procedures ----
;************************************************************************************

;- --- Heightmaps generation ---
;************************************************************************************
; Name: makeDiamondSquare
; Purpose: Creates an array of fractal noise (looks like Perlin) that can be used as
;          a heightmap
; Parameters:
;   - size of the array (a square of size x size)
;   - dispersion
;   - seed: albeit the noise is random, a seed always produces the same result
; Return-value: none, but the result is stored in the "diamond" array
;************************************************************************************
Procedure makeDiamondSquare(gridSize.i,dispersion.f,seed.i)

  Protected gridstep.i,mid.i,n.i
  Protected i.i,i1.i,i2.i
  Protected j.i,j1.i,j2.i,js.i
  Protected u.f,average.f
  Protected dispStep.f,min.f=999999,max.f=-999999,ratio.f
  Protected Dim temp.f(gridSize,gridSize)
  Protected t.i = ElapsedMilliseconds()

  RandomSeed(seed)

  gridstep = gridSize
  dispStep = dispersion

  ; main loop
  While gridstep > 1
    mid=gridstep / 2

    ; Diamond step - calculates new diamond corners from squares
    i=mid
    i1 = i - mid
    i2 = i + mid
    
    While i < gridSize
      j = mid
      j1 = j - mid
      j2 = j + mid
      
      While j < gridSize
        ; Average of surrounding points
        average = (temp(i1,j1)+temp(i1,j2)+temp(i2,j1)+temp(i2,j2))/4
        
        ; calculate random values between -1 and 1
        u = Random(16384)/8192 - 1
        
        ; Diamond value
        temp(i,j) = average + (u * dispStep)
        j + gridstep
        j1 + gridstep
        j2 + gridstep
      Wend
      
      i + gridstep
      i1 + gridstep
      i2 + gridstep
    Wend

    ; square Step - calculates new square corners from diamonds
    i = 0
    i1 = i - mid
    i2 = i + mid
    js = 0
    
    While i < gridSize
      js = mid - js ; toggle start values of j loop
      j = js
      j1 = j - mid
      j2 = j + mid
      
      While j < gridSize ;+1
        average = 0
        If i1 < 0  ; check for need to wrap around i value
          average + temp(i2,j) + temp(i2,j)
        Else
          If i2 > gridSize
            average + temp(i1,j) + temp(i1,j)
          Else
            average + temp(i1,j) + temp(i2,j)
          EndIf
        EndIf
        
        If j1<0  ; check for need to wrap around j value
          average + temp(i,j2) + temp(i,j2)
        Else
          If j2 > gridSize
            average + temp(i,j1) + temp(i,j1)
          Else
            average + temp(i,j1) + temp(i,j2)
          EndIf
        EndIf     
        average/4
        
        ; calculate random value between -1 And 1
        u = Random(16384)/8192 - 1
        
        temp(i,j) = average + (u * dispStep)
        temp(gridSize,j) = temp(0,j) ; copy opposite edge
        j + gridstep
        j1 + gridstep
        j2 + gridstep
        
      Wend
      If j = gridSize
        temp(i,j) = temp(i,0) ; copy opposite edge
      EndIf
      i + mid
      i1 + mid
      i2 + mid
    Wend

    dispStep/2
    gridstep/2
  Wend

; Keep values in Byte range (between 0 and 255)
  For i=0 To gridSize-1
    For j=0 To gridSize-1
      ; Min /max values
      If temp(i,j)<min
        min=temp(i,j)
      Else
        If temp(i,j)>max
          max=temp(i,j)
        EndIf
      EndIf   
    Next j
  Next i
  
  ratio=255.0/(max-min)
  For i=0 To gridSize-1
    For j=0 To gridSize-1
      n=(temp(i,j)-min)*ratio
      If n<0
        diamond(i,j)=0
      Else
        If n>255
          diamond(i,j)=255
        Else
          diamond(i,j)=n
        EndIf
      EndIf
    Next j
  Next i
  
  Debug "Diamond square: "+Str(ElapsedMilliseconds()-t)+" ms."

EndProcedure

;************************************************************************************
; Name: makeVoronoi
; Purpose: Creates a modified Voronoi diagram that can be used to simulate tectonic
;          plates
; Parameters:
;   - size of the array (a square of size x size)
;   - number of regions
;   - number of points in each region (a point is the "center" of a Voronoi plate)
;   - seed: albeit the diagram is random, a seed always produces the same result
; Return-value: none, but the result is stored in the "voronoi" array
;************************************************************************************
Procedure makeVoronoi(gridSize.i,nbregion.i,nbpts.i,seed.i)

  Protected NewList ftpoint.VoronoiPoint2D_struct()
  Protected Dim coef.f(nbregion*nbregion*nbpts)
  Protected i.i,j.i,k.i,n.i
  Protected min.f=999999,max.f=-999999,ratio.f
  Protected Dim temp.f(gridSize,gridSize)
  Protected t.i = ElapsedMilliseconds()
  
  RandomSeed(seed)
  coef(0)=-1
  coef(1)=-2
  coef(2)=0.5
  
  ; In each region of the grid, place some points
  For i=0 To nbregion-1
    For j =0 To nbregion-1
      n=Random(nbpts)
      For k=1 To n
        AddElement(ftpoint())
        ftpoint()\x=i*(gridSize/nbregion)+Random(gridSize/nbregion)
        ftpoint()\y=j*(gridSize/nbregion)+Random(gridSize/nbregion)
        ;Debug StrF(ftpoint()\x)+","+StrF(ftpoint()\y)
      Next k
    Next j
  Next i
  
  ; For each cell of the grid, compute distance to each point, sort points according to distance.
  ; => value = coef * distance
  For i=0 To gridSize-1
    For j=0 To gridSize-1
    
      ForEach ftpoint()     
        ftpoint()\dist = Sqr((ftpoint()\x-i)*(ftpoint()\x-i)+(ftpoint()\y-j)*(ftpoint()\y-j))
      Next
      SortStructuredList(ftpoint(),0,OffsetOf(VoronoiPoint2D_struct\dist),#PB_Sort_Float)
      
      k=0
      ForEach ftpoint()
        temp(i,j) + coef(k)*ftpoint()\dist
        k+1
      Next     
      
      ; Min /max values
      If temp(i,j)<min
        min=temp(i,j)
      Else
        If temp(i,j)>max
          max=temp(i,j)
        EndIf
      EndIf   
      
    Next j
  Next i

  ; Keep values in Byte range (between 0 and 255)
  ratio=255.0/(max-min)
  For i=0 To gridSize-1
    For j=0 To gridSize-1

      n=(temp(i,j)-min)*ratio
      If n<0
        voronoi(i,j)=0
      Else
        If n>255
          voronoi(i,j)=255
        Else
          voronoi(i,j)=n
        EndIf
      EndIf

    Next j
  Next i
  
  Debug "Voronoi: "+Str(ElapsedMilliseconds()-t)+" ms."
EndProcedure

;************************************************************************************
; Name: MakeHeightmap
; Purpose: Mix a diamond-square array and a Voronoi diagram to obtain the final heightmap
; Parameters:
;   - number of the heightmap
;   - size of the array (a square of size x size); Diamond, Voronoi, and Heightmap must
;     have the same size
;   - use the diamond-square array or not
;   - use the Voronoi diagram or not 
; Return-value:  none, but the result is stored in the "heightmap" array
;************************************************************************************
Procedure makeHeightmap(numHmap.i,gridSize.i,Diamond.b,Voronoi.b)

  Protected i.i,j.i,k1.i,k2.i
  Protected min.f=999999,max.f=-999999,ratio.f,temp.f
  Protected Dim temp.f(gridSize,gridSize)
  Protected t.i = ElapsedMilliseconds()
  
  ; for each point of the array:
  For i=0 To gridSize-1
    For j=0 To gridSize-1
    
      ; both arrays will be used
      If Voronoi=#True And Diamond=#True 
        k1=(diamond(i,j) & 255)
        k2=(voronoi(i,j) & 255)
        temp(i,j) = k1*k2
      Else
        ; only diamond-square
        If Diamond=#True
          k1=(diamond(i,j) & 255)
          temp(i,j) = k1
        Else
          ; only Voronoi
          If Voronoi=#True
            k2=(voronoi(i,j) & 255)
            temp(i,j) = k2
          EndIf
        EndIf
      EndIf    
    
      ; Store Min / Max values
      If temp(i,j)<min
        min=temp(i,j)
      Else
        If temp(i,j)>max
          max=temp(i,j)
        EndIf
      EndIf
    
    Next j
  Next i
  
  ; The heightmap values must be between 0 and 255
  ratio = 255.0/(max-min)
  For i=0 To gridSize-1
    For j=0 To gridSize-1
      temp = (temp(i,j)-min)*ratio
      ;Debug temp
      If temp<0
        heightmap(numHmap,i,j)=0
      Else
        If temp>255
          heightmap(numHmap,i,j)=255
        Else
          heightmap(numHmap,i,j)=temp
        EndIf
      EndIf
    Next j
  Next i
  
  Debug "Make heightmap: "+Str(ElapsedMilliseconds()-t)+" ms."

EndProcedure

;************************************************************************************
; Name: eraseHeightmapSeam
; Purpose: Blends the seam between two heightmaps
; Parameters:
;   - number of the first heightmap
;   - side of the first heightmap to blend (top, bottom, left, or right)
;   - number of the second heightmap
;   - side of the second heightmap to blend (top, bottom, left, or right)
;   - size of the array (a square of size x size); both Heightmaps must
;     have the same size
;   - percentage of that size that will be blended
; Return-value:  none
;************************************************************************************
Procedure eraseHeightmapSeam(numHmap1.i,side1.b,numHmap2.i,side2.b,gridSize.i,percent.b)

  Protected i.i,j.i,width.i,seamwidth.i,color.i
  Protected x1.i,y1.i,x2.i,y2.i
  Protected h1.i,h2.i,h.i
  Protected ratio1.f,ratio2.f
  Protected t.i = ElapsedMilliseconds()
  
  width = gridSize-1
  seamwidth = ( width*percent/100 )/2
    
  For i=0 To seamwidth
  
    ratio1 = (i/2)/seamwidth
    ratio2 = 1.0 - ratio1

    For j=0 To gridSize-1
    
      Select side1
        Case #HEIGHTMAP_TOP
          x1 = j
          y1 = seamwidth-i
        Case #HEIGHTMAP_BOTTOM
          x1 = j
          y1 = width-seamwidth+i
        Case #HEIGHTMAP_LEFT
          x1 = seamwidth-i
          y1 = j
        Case #HEIGHTMAP_RIGHT
          x1 = width-seamwidth+i
          y1 = j
      EndSelect
      
      Select side2
        Case #HEIGHTMAP_TOP
          x2 = j
          y2 = seamwidth-i
          If side1 = #HEIGHTMAP_TOP Or side1 = #HEIGHTMAP_RIGHT
            x2=gridSize-1-j
          EndIf         
        Case #HEIGHTMAP_BOTTOM
          x2 = j
          y2 = width-seamwidth+i
          If side1 = #HEIGHTMAP_BOTTOM Or side1 = #HEIGHTMAP_LEFT
            x2=gridSize-1-j
          EndIf
        Case #HEIGHTMAP_LEFT
          x2 = seamwidth-i
          y2 = j
          If side1 = #HEIGHTMAP_BOTTOM Or side1 = #HEIGHTMAP_LEFT
            x2=gridSize-1-j
          EndIf
        Case #HEIGHTMAP_RIGHT
          x2 = width-seamwidth+i
          y2 = j
          If side1 = #HEIGHTMAP_TOP Or side1 = #HEIGHTMAP_RIGHT
            x2=gridSize-1-j
          EndIf         
      EndSelect
    
      h1=heightmap(numHmap1,x1,y1)
      h2=heightmap(numHmap2,x2,y2)

      h=h1*ratio2 + h2*ratio1      
      heightmap(numHmap1,x1,y1)=h
      
      h=h2*ratio2 + h1*ratio1
      heightmap(numHmap2,x2,y2)=h      
    Next j
    
  Next i
  
  Debug "Erase seams: "+Str(ElapsedMilliseconds()-t)+" ms."

EndProcedure


;************************************************************************************
; Name: ComputeGeoStats
; Purpose: Compute height, slope, temperature, ice/snow, terrain type and color for
;          each cell of an heightmap
; Parameters:
;   - number of the heightmap
;   - size of the array (a square of size x size)
;   - sea level (between 0-255)
;   - Temperature at poles (usually cold)
;   - Temperature at equator (usually hot)
;   - seed: albeit the map is random, a seed always produces the same result
; Return-value:  none, but the results are stored in the "geostat" array
;************************************************************************************
Procedure computeGeoStats(numHmap.i,gridSize.i,waterLevel.i,tpol.f,teq.f,seed.i)

  Protected heightAboveWater.i
  Protected equivHeightToTemp.f = -0.4; temperature loss when you climb up 1 unit
  Protected i.i,j.i,k.i
  Protected halfGridSize.i = gridSize/2
  Protected x.i,y.i,newx.i,newy.i
  Protected p.i,h1.i,h2.i,maxp.i
  Protected red1.i,red2.i,green1.i,green2.i,blue1.i,blue2.i,red.f,green.f,blue.f
  Protected exp.f,cumul.f,icegradient.f,icecolor.i
  Protected vector3.Vector3_struct
  Protected t.i = ElapsedMilliseconds()

  RandomSeed(seed)
  
  gridSize-1
  For i=0 To gridSize
    For j = 0 To gridSize
    
      ; Latitude
      If numHmap=#TOPSIDE Or numHmap=#BOTTOMSIDE
        vector3\x = i - halfGridSize
        vector3\y = halfGridSize
        vector3\z = j - halfGridSize
    
        VEC3_NORMALIZE(vector3)
        
        vector3\x * halfGridSize
        vector3\y * halfGridSize
        vector3\z * halfGridSize

        geostat(numHmap,i,j)\latitude = Cos(vector3\y/gridSize * #PI)
      Else
        vector3\x = i - halfGridSize
        vector3\y = j - halfGridSize
        vector3\z = halfGridSize
    
        VEC3_NORMALIZE(vector3)
        
        vector3\x * halfGridSize
        vector3\y * halfGridSize
        vector3\z * halfGridSize

        geostat(numHmap,i,j)\latitude = Cos(vector3\y/gridSize * #PI)
      EndIf
    
      ; Altitude
      geostat(numHmap,i,j)\height = heightmap(numHmap,i,j)
      
      ; Slope
      h1= geostat(numHmap,i,j)\height
      geostat(numHmap,i,j)\slope=0
      
      If i>0
        h2=heightmap(numHmap,i-1,j)
        p=Abs( h1 - h2 )
        If p>geostat(numHmap,i,j)\slope
          geostat(numHmap,i,j)\slope = p
        EndIf
      EndIf

      If i<gridSize
        h2=heightmap(numHmap,i+1,j)
        p=Abs( h1 - h2 )
        If p>geostat(numHmap,i,j)\slope
          geostat(numHmap,i,j)\slope = p
        EndIf
      EndIf

      If j>0
        h2=heightmap(numHmap,i,j-1)
        p=Abs( h1 - h2 )
        If p>geostat(numHmap,i,j)\slope
          geostat(numHmap,i,j)\slope = p
        EndIf
      EndIf

      If j<gridSize
        h2=heightmap(numHmap,i,j+1)
        p=Abs( h1 - h2 )
        If p>geostat(numHmap,i,j)\slope
          geostat(numHmap,i,j)\slope = p
        EndIf
      EndIf
      If geostat(numHmap,i,j)\slope > maxp
        maxp = geostat(numHmap,i,j)\slope
      EndIf

      ; Temperature
      heightAboveWater=geostat(numHmap,i,j)\height - waterlevel
      If heightAboveWater <0
        heightAboveWater = 0
      EndIf
      geostat(numHmap,i,j)\temperature = tpol + geostat(numHmap,i,j)\latitude*(teq-tpol) + equivHeightToTemp*heightAboveWater
      
      If geostat(numHmap,i,j)\height<waterlevel
        ; This point is under sea level
        geostat(numHmap,i,j)\isWater=2
      
        ; underwater points => shades of blue
        red1=0:green1=3:blue1=92
        red2=0:green2=130:blue2=220
        
        red=red1+(red2-red1)*(geostat(numHmap,i,j)\height/waterlevel)
        green=green1+(green2-green1)*(geostat(numHmap,i,j)\height/waterlevel)
        blue=blue1+(blue2-blue1)*(geostat(numHmap,i,j)\height/waterlevel)
      
        ; water depth causes temperature variation
        geostat(numHmap,i,j)\temperature-(geostat(numHmap,i,j)\height-waterlevel)*0.075
      Else       
        ; This point is above sea level
        geostat(numHmap,i,j)\isWater=0
        
        ; Probability computation for each terrain type
        cumul=0:red=0:green=0:blue=0
        For k=1 To 5
          exp = Abs(terrain(k)\temperatureAv - geostat(numHmap,i,j)\temperature)/terrain(k)\temperatureVt
          exp + Abs(terrain(k)\slopeAv - geostat(numHmap,i,j)\slope)/terrain(k)\slopeVt
          exp + Abs(terrain(k)\heightAv*waterLevel - geostat(numHmap,i,j)\height)/terrain(k)\heightVt
          geostat(numHmap,i,j)\probaterrain[k] = Exp(-exp); Pow(#E,-exp)
          
          cumul + geostat(numHmap,i,j) \ probaterrain[k]
        Next k
        
        ; Mix terrains colors according to probability
        For k=1 To 5
          geostat(numHmap,i,j)\probaterrain[k] / cumul
          blue+Red(terrain(k)\colour)*geostat(numHmap,i,j)\probaterrain[k]
          green+Green(terrain(k)\colour)*geostat(numHmap,i,j)\probaterrain[k]
          red+Blue(terrain(k)\colour)*geostat(numHmap,i,j)\probaterrain[k]
        Next k
                
      EndIf
  
      geostat(numHmap,i,j)\colour = RGB(red,green,blue)
      
    Next j
  Next i
  
  ; Generate a few rivers
  j=0
  For i=1 To Random(6)
    Repeat
      x=Random(gridSize)
      y=Random(gridSize)
      j+1
    Until geostat(numHmap,x,y)\isWater<>2 Or j=1000
    
    Repeat
      geostat(numHmap,x,y)\colour = RGB(0,130,220)
      geostat(numHmap,x,y)\isWater=1
      h1=999
      newx=0:newy=0
      If x>0
        If geostat(numHmap,x-1,y)\isWater <>1 And geostat(numHmap,x-1,y)\height < h1
          newx=x-1:newy=y
          h1=geostat(numHmap,newx,newy)\height
        EndIf 
      EndIf
      If x<gridSize
        If geostat(numHmap,x+1,y)\isWater <>1 And geostat(numHmap,x+1,y)\height < h1
          newx=x+1:newy=y
          h1=geostat(numHmap,newx,newy)\height
        EndIf 
      EndIf
      If y>0
        If geostat(numHmap,x,y-1)\isWater <>1 And geostat(numHmap,x,y-1)\height < h1
          newx=x:newy=y-1
          h1=geostat(numHmap,newx,newy)\height
        EndIf 
      EndIf
      If y<gridSize
        If geostat(numHmap,x,y+1)\isWater <>1 And geostat(numHmap,x,y+1)\height < h1
          newx=x:newy=y+1
          h1=geostat(numHmap,newx,newy)\height
        EndIf 
      EndIf
      x=newx:y=newy
      
    Until geostat(numHmap,x,y)\isWater=2 Or (newx=0 And newy=0)
    
  Next i

  ; Add ice/snow, slope marking
  For i=0 To gridSize
    For j=0 To gridSize
      red=Red(geostat(numHmap,i,j)\colour)
      green=Green(geostat(numHmap,i,j)\colour)
      blue=Blue(geostat(numHmap,i,j)\colour)
    
      ; Ice / snow
      If geostat(numHmap,i,j)\temperature < 0
        If geostat(numHmap,i,j)\temperature < 0
          icegradient= -geostat(numHmap,i,j)\temperature / 5
        EndIf
        
        If icegradient>1
          icegradient=1
        EndIf         
        icecolor=247+Random(8)
        
        red=red*(1-icegradient) + icecolor*icegradient
        green=green*(1-icegradient) + icecolor*icegradient
        blue=blue*(1-icegradient) + icecolor*icegradient        
      EndIf 
      
      ; Slopes are darker
      If geostat(numHmap,i,j)\isWater<2    
        red-geostat(numHmap,i,j)\slope * 3
        green-geostat(numHmap,i,j)\slope * 3
        blue-geostat(numHmap,i,j)\slope * 3        
      EndIf     
        
      ; RGB must stay between 0-255
      If red<0 
        red=0
      Else
        If red>255
          red=255
        EndIf
      EndIf
      If green<0 
        green=0
      Else
        If green>255
          green=255
        EndIf
      EndIf
      If blue<0 
        blue=0
      Else
        If blue>255
          blue=255
        EndIf
      EndIf
        
      geostat(numHmap,i,j)\colour = RGB(blue,green,red)

    Next j
  Next i
  
  Debug "Compute geostats: "+Str(ElapsedMilliseconds()-t)+" ms."

EndProcedure

;************************************************************************************
; Name: MakeTerrainImage
; Purpose: Draw a map from heightmap and geostats arrays
; Parameters:
;   - number of the heightmap / geostats array
;   - size of the array (a square of size x size)
;   - Type of map to draw: "terrain","altitude","slope", "temperature"
; Return-value:  generated image number
;************************************************************************************
Procedure.i makeTerrainImage(numHmap.i,gridSize.i,imageType.b)
  
  Protected numImage.i
  Protected i.i,j.i,color.i,icecolor.i
  Protected red.f,green.f,blue.f,icegradient.f
  Protected t.i =ElapsedMilliseconds()
  
  numImage = CreateImage(#PB_Any,gridSize,gridSize)
  
  CompilerIf #USEAPIDRAWING = #True
    Protected hDC.i,hBmp.i,*mem.i,picl_X.i,picl_Y.i,picl_D.i
    Protected bmi.myBITMAPINFO,*pixel.LONG
  
    hDC=StartDrawing(ImageOutput(numImage))
    hBmp = ImageID(numImage)
    picl_X = ImageWidth(numImage)
    picl_Y = ImageHeight(numImage)+1
    picl_D = ImageDepth(numImage)  
    
    *mem = AllocateMemory(picl_X*picl_Y*4)
    bmi.myBITMAPINFO
    bmi\bmiHeader\biSize   = SizeOf(BITMAPINFOHEADER)
    bmi\bmiHeader\biWidth  = picl_X
    bmi\bmiHeader\biHeight = picl_Y
    bmi\bmiHeader\biPlanes = 1
    bmi\bmiHeader\biBitCount = 32
    bmi\bmiHeader\biCompression = #BI_RGB
    GetDIBits_(hDC,hBmp,1,picl_Y,*mem,bmi,#DIB_RGB_COLORS)
  
    *pixel = *mem
  CompilerElse
    StartDrawing(ImageOutput(numImage))
  CompilerEndIf
  
  For j=gridSize-1 To 0 Step-1
    For i=0 To gridSize-1
    
      Select imageType
      
        Case #MAPTYPE_ALTITUDE
          red = geostat(numHmap,i,j)\height
          green = red:blue = Red ; (greyscale)

          
        Case #MAPTYPE_TEMPERATURE
          If i>1 And i<gridSize-1 And j>1 And j<gridSize-1 And geostat(numHmap,i,j)\isWater>0 And (geostat(numHmap,i-1,j)\isWater=0 Or geostat(numHmap,i+1,j)\isWater=0 Or geostat(numHmap,i,j-1)\isWater=0 Or geostat(numHmap,i,j+1)\isWater=0)
            color=0
          Else
            red = (geostat(numHmap,i,j)\temperature - tpol)*255/(teq-tpol)
            If red>255
              red=255
            EndIf
            If red<0
              red=0
            EndIf
            blue = 255-red
            green=blue/2
          EndIf
          
        Case #MAPTYPE_SLOPE
          red = geostat(numHmap,i,j)\slope*32
          If red > 255
            red = 255
          EndIf
          green = red:blue = Red ; (greyscale)
          
        Case #MAPTYPE_TERRAIN
          ; Add ice/snow, slope marking
            red=Red(geostat(numHmap,i,j)\colour)
            green=Green(geostat(numHmap,i,j)\colour)
            blue=Blue(geostat(numHmap,i,j)\colour)
          
            ; Ice / snow
            If geostat(numHmap,i,j)\temperature < 0
              If geostat(numHmap,i,j)\temperature < 0
                icegradient= -geostat(numHmap,i,j)\temperature / 5
              EndIf
              
              If icegradient>1
                icegradient=1
              EndIf         
              icecolor=247+Random(8)
              
              red=red*(1-icegradient) + icecolor*icegradient
              green=green*(1-icegradient) + icecolor*icegradient
              blue=blue*(1-icegradient) + icecolor*icegradient        
            EndIf 
            
            ; Slopes are darker
            If geostat(numHmap,i,j)\isWater<2    
              red-geostat(numHmap,i,j)\slope * 4
              green-geostat(numHmap,i,j)\slope * 4
              blue-geostat(numHmap,i,j)\slope * 4        
            EndIf     
              
            ; RGB must stay between 0-255
            If red<0 
              red=0
            Else
              If red>255
                red=255
              EndIf
            EndIf
            If green<0 
              green=0
            Else
              If green>255
                green=255
              EndIf
            EndIf
            If blue<0 
              blue=0
            Else
              If blue>255
                blue=255
              EndIf
            EndIf
          
      EndSelect
      
      CompilerIf #USEAPIDRAWING = #True
        color=RGB(red,green,blue)
        *pixel\l=color
        *pixel + 4
      CompilerElse
        color=RGB(blue,green,red)
        Plot(i,j,color)
      CompilerEndIf
    Next i
  Next j
  
  CompilerIf #USEAPIDRAWING = #True
    SetDIBits_(hDC,hBmp,1,gridSize+1,*mem,bmi,#DIB_RGB_COLORS)
  CompilerEndIf
  
  StopDrawing()
  
  CompilerIf #USEAPIDRAWING = #True
    FreeMemory(*mem)
  CompilerEndIf
  
  Debug "Make terrain image: "+Str(ElapsedMilliseconds()-t)+" ms."

  ProcedureReturn numImage
EndProcedure



;- --- Mesh generation ---
;************************************************************************************
; Name: initializeMeshDatas
; Purpose: Creates a list of vertices and a list of triangles for a "plane" mesh
; Parameters:
;   - vertices list
;   - triangles list
;   - sizes x,y,z of the plane
; Return-value: none, but the results are stored in the "vertex" and "triangle" lists()
;************************************************************************************
Procedure.i initializeMeshDatas(List vertex.vertex_struct(),List triangle.triangle_struct(),sizeX.f=1,sizeY.f=1,sizeZ.f=1)

  Protected nbVert.i , nbTri.i              ; Number of vertices and faces
  Protected x.f,y.f,z.f                     ; vertex position
  Protected i.i
  
  ; Read number of vertices and triangles
  Read.i nbVert
  Read.i nbTri
  
  ; Read and store vertices position, normals, uv coords
  For i = 1 To nbVert
    AddElement(vertex())
   
    Read.f x
    Read.f y
    Read.f z
   
    vertex()\x = x * sizex
    vertex()\y = y * sizey
    vertex()\z = z * sizez
   
    Read.f vertex()\nx
    Read.f vertex()\ny
    Read.f vertex()\nz
    Read.f vertex()\u
    Read.f vertex()\v
   
    vertex()\colour = $FFFFFF
    
  Next i   
 
  ;Read and store faces infos
  For i=1 To nbTri
    AddElement(triangle())
   
    Read.i triangle()\v1
    Read.i triangle()\v2
    Read.i triangle()\v3
   
    ; store vertices' pointers, for faster access
    SelectElement(vertex(),triangle()\v1)
    triangle()\ptrV1 = @vertex()
    SelectElement(vertex(),triangle()\v2)
    triangle()\ptrV2 = @vertex()
    SelectElement(vertex(),triangle()\v3)
    triangle()\ptrV3 = @vertex()
  Next i
EndProcedure

;************************************************************************************
; Name: CopyAndRotateVertexList
; Purpose: Copy and rotate a vertex list to make all the sides of the planet from
;          the "top" one. (faster than tesselate the six sides)
; Parameters:
;   - vertices list for "top" side
;   - vertices list to store the result
;   - side you want to create
; Return-value: none, but the result is stored in the "vertex" list()
;************************************************************************************
Procedure copyAndRotateVertexList(List vertex.vertex_struct(),List newVertex.vertex_struct(),side.i)  
  ; Note: This procedure assumes the "vertex" list is the "TOP" one
  CopyList(vertex(),newVertex())
  Select side
    Case #FRONTSIDE
      ForEach newVertex()
        Swap newVertex()\y,newVertex()\z
        newVertex()\y = -newVertex()\y        
      Next newVertex()

    Case #BACKSIDE
      ForEach newVertex()
        Swap newVertex()\y,newVertex()\z
        newVertex()\y = -newVertex()\y
        
        newVertex()\x = -newVertex()\x
        newVertex()\z = -newVertex()\z
      Next newVertex()

    Case #LEFTSIDE
      ForEach newVertex()
        Swap newVertex()\y,newVertex()\z
        newVertex()\y = -newVertex()\y        
        
        Swap newVertex()\x,newVertex()\z
        newVertex()\x = -newVertex()\x
      Next newVertex()

    Case #RIGHTSIDE
      ForEach newVertex()
        Swap newVertex()\y,newVertex()\z
        newVertex()\y = -newVertex()\y        
        
        Swap newVertex()\x,newVertex()\z
        newVertex()\z = -newVertex()\z
      Next newVertex()

    Case #BOTTOMSIDE
      ForEach newVertex()
        newVertex()\y = -newVertex()\y        
        newVertex()\z = -newVertex()\z        
      Next newVertex()
      
  EndSelect
  
EndProcedure

;************************************************************************************
; Name: tesselateMesh
; Purpose: Break each triangle of the mesh into four smaller triangle to increase LOD
; Parameters:
;   - vertices list of the mesh
;   - triangles list of the mesh
;   - number of times you want to repeat the operation
; Return-value: none, but the results are stored in the "vertex" and "triangle" lists()
;************************************************************************************
Procedure.i tesselateMesh(List vertex.vertex_struct(),List triangle.triangle_struct(),depth.i)

  Protected i.i,j.i
  Protected t0.i
 
  Protected NewList newVertex.vertex_struct()
  Protected NewList newTriangle.triangle_struct()
 
  Protected v1_2.vertex_struct,v2_3.vertex_struct,v1_3.vertex_struct
  Protected vertNum1.i,vertNum2.i,vertNum3.i,vertNum1_2.i,vertNum1_3.i,vertNum2_3.i
  
  Protected NewMap antidup.i() ; ugly fix to avoid duplicate vertices...
  
  
  ; Tesselate!
  t0 = ElapsedMilliseconds()
  For i=1 To depth
    ClearList(newVertex())
    ClearList(newTriangle())
    ClearMap(antidup())
   
    ForEach triangle()
      ; Compute new points (in the middle of each side of the triangle)
      MIDDLEPOINT(v1_2,triangle()\ptrV1,triangle()\ptrV2)
      MIDDLEPOINT(v1_3,triangle()\ptrV1,triangle()\ptrV3)
      MIDDLEPOINT(v2_3,triangle()\ptrV2,triangle()\ptrV3)
     
      ; store them in the "new vertex" list, and check if the already exists
      CREATENEWVERTEX(triangle()\ptrV1,vertNum1)
      CREATENEWVERTEX(triangle()\ptrV2,vertNum2)
      CREATENEWVERTEX(triangle()\ptrV3,vertNum3)
      CREATENEWVERTEX(v1_2,vertNum1_2)
      CREATENEWVERTEX(v1_3,vertNum1_3)
      CREATENEWVERTEX(v2_3,vertNum2_3)
     
      ; Create the new triangles
      AddElement(newTriangle())
      newTriangle()\v1 = vertNum1
      newTriangle()\v2 = vertNum1_2
      newTriangle()\v3 = vertNum1_3

      AddElement(newTriangle())
      newTriangle()\v1 = vertNum1_2
      newTriangle()\v2 = vertNum2
      newTriangle()\v3 = vertNum2_3
     
      AddElement(newTriangle())
      newTriangle()\v1 = vertNum1_3
      newTriangle()\v2 = vertNum2_3
      newTriangle()\v3 = vertNum3

      AddElement(newTriangle())
      newTriangle()\v1 = vertNum2_3
      newTriangle()\v2 = vertNum1_3
      newTriangle()\v3 = vertNum1_2
    Next triangle()
   
    ; The new list replaces the old one
    ClearList(vertex())
    ClearList(triangle())
   
    CopyList(newVertex(),vertex())
    CopyList(newTriangle(),triangle())
   
    ; Store vertices' pointers, for faster access
    ForEach triangle()
      SelectElement(vertex(),triangle()\v1)
      triangle()\ptrV1 = @vertex()
      SelectElement(vertex(),triangle()\v2)
      triangle()\ptrV2 = @vertex()
      SelectElement(vertex(),triangle()\v3)
      triangle()\ptrV3 = @vertex()      
    Next triangle()   
    
  Next i 
  
  FreeMap(antidup())
  
  Debug "Depth = " + Str(depth)
  Debug "Nb vertices = " + Str(ListSize(vertex()))
  Debug "Nb triangles = " + Str(ListSize(triangle()))
  Debug "Mesh tesselation: "+Str(ElapsedMilliseconds()-t0) + " ms."
  
EndProcedure


;************************************************************************************
; Name: applyHeightmap
; Purpose: Apply a heightmap to a side of the planet
; Parameters:
;   - Heightmap number
;   - vertices list of the mesh
;   - planet size
; Return-value: none, but the result is stored in the "vertex" list()
;************************************************************************************
Procedure applyHeightmap(numHmap.i,List vertex.vertex_struct(),planetRadius.f,reliefExageration.f = 0.2,gridSize.i = #GRIDSIZE,waterLevel.i=#WATERLEVEL)
  Protected height.i
  Protected displacement.f
  
  gridSize-1
  ForEach vertex()
    height = heightmap(numHmap,Int(vertex()\u * gridSize),Int(vertex()\v * gridSize)) - waterLevel
    If height > 0
      displacement = height*reliefExageration
    Else
      displacement = 0
    EndIf
    
    ; This transforms the cube into a sphere
    VEC3_NORMALIZETOLENGTH(vertex(),planetRadius+displacement)
    
  Next vertex()
  
EndProcedure

;************************************************************************************
; Name: weldSideBorders
; Purpose: Adjust the vertices' coords on the borders of the planet "sides", to avoid
;          gaps between the submeshes
; Parameters:
;   none
; Return-value: none, but the results are stored in the "vertex" and "triangle" lists()
;************************************************************************************
Procedure weldSideBorders()
  Protected i.i
  
  ; First , Find and store the vertices which are along the borders of the mesh
  For i=0 To 5
    ForEach planetSide(i)\vertex()      
      If planetSide(i)\vertex()\v = 0
        AddElement(planetSide(i)\topBorder())
        planetSide(i)\topBorder()\ptrVert = @planetSide(i)\vertex()
        planetSide(i)\topBorder()\u = planetSide(i)\vertex()\u
        planetSide(i)\topBorder()\v = planetSide(i)\vertex()\v
      EndIf
      If planetSide(i)\vertex()\v = 1
        AddElement(planetSide(i)\bottomBorder())
        planetSide(i)\bottomBorder()\ptrVert = @planetSide(i)\vertex()
        planetSide(i)\bottomBorder()\u = planetSide(i)\vertex()\u
        planetSide(i)\bottomBorder()\v = planetSide(i)\vertex()\v
      EndIf
      If planetSide(i)\vertex()\u = 0
        AddElement(planetSide(i)\leftBorder())
        planetSide(i)\leftBorder()\ptrVert = @planetSide(i)\vertex()
        planetSide(i)\leftBorder()\u = planetSide(i)\vertex()\u
        planetSide(i)\leftBorder()\v = planetSide(i)\vertex()\v
      EndIf
      If planetSide(i)\vertex()\u = 1
        AddElement(planetSide(i)\rightBorder())
        planetSide(i)\rightBorder()\ptrVert = @planetSide(i)\vertex()
        planetSide(i)\rightBorder()\u = planetSide(i)\vertex()\u
        planetSide(i)\rightBorder()\v = planetSide(i)\vertex()\v
      EndIf
    Next planetSide(i)\vertex()
  Next i
   
  ; Weld top to front/back/left/right
  WELDBORDERVERTICES(#TOPSIDE,bottom,u,#PB_Sort_Ascending,#FRONTSIDE,top,u)
  WELDBORDERVERTICES(#TOPSIDE,left,v,#PB_Sort_Ascending,#LEFTSIDE,top,u)
  WELDBORDERVERTICES(#TOPSIDE,right,v,#PB_Sort_Descending,#RIGHTSIDE,top,u)
  WELDBORDERVERTICES(#TOPSIDE,top,u,#PB_Sort_Descending,#BACKSIDE,top,u)
  
  ; Weld bottom to front/back/left/right
  WELDBORDERVERTICES(#BOTTOMSIDE,top,u,#PB_Sort_Ascending,#FRONTSIDE,bottom,u)
  WELDBORDERVERTICES(#BOTTOMSIDE,left,v,#PB_Sort_Descending,#LEFTSIDE,bottom,u)
  WELDBORDERVERTICES(#BOTTOMSIDE,right,v,#PB_Sort_Ascending,#RIGHTSIDE,bottom,u)
  WELDBORDERVERTICES(#BOTTOMSIDE,bottom,u,#PB_Sort_Descending,#BACKSIDE,bottom,u)
  
  ; Weld front/back/left/right
  WELDBORDERVERTICES(#FRONTSIDE,right,v,#PB_Sort_Ascending,#RIGHTSIDE,left,v)
  WELDBORDERVERTICES(#RIGHTSIDE,right,v,#PB_Sort_Ascending,#BACKSIDE,left,v)
  WELDBORDERVERTICES(#BACKSIDE,right,v,#PB_Sort_Ascending,#LEFTSIDE,left,v)
  WELDBORDERVERTICES(#LEFTSIDE,right,v,#PB_Sort_Ascending,#FRONTSIDE,left,v)
  
  ; Clear the lists to save memory
  For i=0 To 5
    ClearList(planetSide(i)\topBorder())
    ClearList(planetSide(i)\bottomBorder())
    ClearList(planetSide(i)\leftBorder())
    ClearList(planetSide(i)\rightBorder())
  Next i
  
  
EndProcedure


;************************************************************************************
; Name: createMeshFromLists
; Purpose: Create a mesh from the six "vertex()" and "triangle()" lists of the planet.
;          (There's one submesh for each side of the planet)
; Parameters:
;   none
; Return-value: Planet mesh number
;************************************************************************************
Procedure.i createMeshFromLists(Array planetSide.planetSide_struct(1))
  Protected numMesh.i,i.i
  
  ; Create the mesh
  numMesh = CreateMesh(#PB_Any)
  
  ; Add vertices to mesh
  For i =0 To 5
    AddSubMesh()
    
    ForEach planetSide(i)\vertex()   
      AddMeshVertex(planetSide(i)\vertex()\x,planetSide(i)\vertex()\y,planetSide(i)\vertex()\z)
      ;MeshVertexColor(vertex()\colour)
      ;MeshVertexNormal(vertex()\nx,vertex()\ny,vertex()\nz)
      MeshVertexTextureCoordinate(planetSide(i)\vertex()\u,planetSide(i)\vertex()\v)
    Next planetSide(i)\vertex()
   
    ; Add triangles to mesh
    ForEach planetSide(i)\triangle()
      AddMeshFace(planetSide(i)\triangle()\v1,planetSide(i)\triangle()\v2,planetSide(i)\triangle()\v3)
    Next planetSide(i)\triangle()        
    
  Next i
  
  ; Finalize mesh
  FinishMesh()
  NormalizeMesh(numMesh) ; NB: it doesn't normalize between submeshes => visible seams
  
  ; Default material = planet sides textures 
  For i =0 To 5
    SetMeshMaterial(numMesh, MaterialID(planetSide(i)\numMaterial),i)
  Next i
  
  ProcedureReturn numMesh
 
EndProcedure   

;************************************************************************************
; Name: createPlanet
; Purpose: Create a planet entity from the "planetSide" array, which contains heighmaps
;          and textures for the six "sides" of the planet  , then attach it to a node    
; Parameters:
;   - "PlanetSide" array
;   - Parent node number
;   - Tesselation level
; Return-value: Planet mesh number
;************************************************************************************
Procedure createPlanet(Array planetSide.planetSide_struct(1), tesselationLevel.i)
  Protected i.i
  Protected planetMesh.i,planetEntity.i
  
  
 ; Clear vertices and triangles list
  For i=0 To 5
    ClearList(planetSide(i)\vertex())
    ClearList(planetSide(i)\triangle())
  Next i

  
  ; Read datas to create the "top" side and tesselate it
  Restore meshDatas
  InitializeMeshDatas(planetSide(#TOPSIDE)\vertex(),planetSide(#TOPSIDE)\triangle(),500,500,500)
  tesselateMesh(planetSide(#TOPSIDE)\vertex(),planetSide(#TOPSIDE)\triangle(),tesselationLevel)
  
  ; Copy it to the other sides
  For i=1 To 5
    CopyAndRotateVertexList(planetSide(#TOPSIDE)\vertex(),planetSide(i)\vertex(),i)
    CopyList(planetSide(#TOPSIDE)\triangle(),planetSide(i)\triangle())
  Next i
  
  ; Apply heightmaps
  For i=0 To 5
    applyHeightmap(i,planetSide(i)\vertex(),500,0.4,#GRIDSIZE,#WATERLEVEL)
  Next i
  
  ; Ensures each side is correctly welded to the next ones
  weldSideBorders()
  
  ; Create the mesh and the entity
  planetMesh = CreateMeshFromLists(planetSide())
  planetEntity = CreateEntity(#PB_Any,MeshID(planetMesh),#PB_Material_None)
  FreeMesh(planetMesh)
  
  ProcedureReturn planetEntity  
EndProcedure 


DisableExplicit
;************************************************************************************

;- ---------- Main program ----------
;- Screen initialisation 
If InitEngine3D() = 0 
   MessageRequester( "Error" , "Can't initialize 3D, check if engine3D.dll is available" , 0 ) 
   End 
EndIf
If InitSprite() = 0 Or InitSprite3D() = 0 Or InitKeyboard() = 0  Or InitMouse() = 0 
   MessageRequester( "Error" , "Can't find DirectX 7.0 or above" , 0 ) 
   End 
EndIf 

OpenWindow(0,0, 0, 800 , 600 ,"Telluric Planet Generator v3",#PB_Window_ScreenCentered) 
OpenWindowedScreen(WindowID(0),0,0, 800 , 600,0,0,0,#PB_Screen_SmartSynchronization) 
KeyboardMode(#PB_Keyboard_International)



;- Heightmaps
; Prepare six heightmaps (one for each side)
dispersion=#GRIDSIZE
t0 = ElapsedMilliseconds()
For i=0 To 5
  ; Compute a "Perlin-like" heightmap
  makeDiamondSquare(#GRIDSIZE,dispersion,#PLANETSEED+i)
  
  ; Voronoi can add a nice "fractured" touch to heightmaps, but can be long to generate
  makeVoronoi(#GRIDSIZE,3,3,#PLANETSEED+i)
  
  ; Pass #true as last parameter to use Voronoi
  makeHeightmap(i,#GRIDSIZE,#True,#False)  
Next i
Debug "Heightmaps generation: "+Str(ElapsedMilliseconds()-t0)+" ms."

;- Blend seams between heightmaps
eraseHeightmapSeam(#TOPSIDE,#HEIGHTMAP_BOTTOM,#FRONTSIDE,#HEIGHTMAP_TOP,#GRIDSIZE,25)
eraseHeightmapSeam(#TOPSIDE,#HEIGHTMAP_TOP,#BACKSIDE,#HEIGHTMAP_TOP,#GRIDSIZE,25)
eraseHeightmapSeam(#TOPSIDE,#HEIGHTMAP_LEFT,#LEFTSIDE,#HEIGHTMAP_TOP,#GRIDSIZE,25) 
eraseHeightmapSeam(#TOPSIDE,#HEIGHTMAP_RIGHT,#RIGHTSIDE,#HEIGHTMAP_TOP,#GRIDSIZE,25)

eraseHeightmapSeam(#FRONTSIDE,#HEIGHTMAP_LEFT,#LEFTSIDE,#HEIGHTMAP_RIGHT,#GRIDSIZE,25)
eraseHeightmapSeam(#FRONTSIDE,#HEIGHTMAP_RIGHT,#RIGHTSIDE,#HEIGHTMAP_LEFT,#GRIDSIZE,25)
eraseHeightmapSeam(#BACKSIDE,#HEIGHTMAP_RIGHT,#LEFTSIDE,#HEIGHTMAP_LEFT,#GRIDSIZE,25)
eraseHeightmapSeam(#BACKSIDE,#HEIGHTMAP_LEFT,#RIGHTSIDE,#HEIGHTMAP_RIGHT,#GRIDSIZE,25)
 
eraseHeightmapSeam(#BOTTOMSIDE,#HEIGHTMAP_TOP,#FRONTSIDE,#HEIGHTMAP_BOTTOM,#GRIDSIZE,25)
eraseHeightmapSeam(#BOTTOMSIDE,#HEIGHTMAP_BOTTOM,#BACKSIDE,#HEIGHTMAP_BOTTOM,#GRIDSIZE,25)
eraseHeightmapSeam(#BOTTOMSIDE,#HEIGHTMAP_LEFT,#LEFTSIDE,#HEIGHTMAP_BOTTOM,#GRIDSIZE,25)
eraseHeightmapSeam(#BOTTOMSIDE,#HEIGHTMAP_RIGHT,#RIGHTSIDE,#HEIGHTMAP_BOTTOM,#GRIDSIZE,25)


; Read terrain types
i=0
Read.s name
While name <> "XXX END XXX"
  i+1
  terrain(i)\name = name
  Read.f terrain(i)\heightAv
  Read.f terrain(i)\heightVt
  Read.f terrain(i)\slopeAv
  Read.f terrain(i)\slopeVt
  Read.f terrain(i)\temperatureAv
  Read.f terrain(i)\temperatureVt
  Read.i terrain(i)\colour
  
  Read.s name
Wend

;- Textures and materials
Add3DArchive(".",#PB_3DArchive_FileSystem)
For i=0 To 5
   computeGeoStats(i,#GRIDSIZE,#WATERLEVEL,tpol,teq,#PLANETSEED)
   
   planetSide(i)\numTexture = makeTerrainImage(i,#GRIDSIZE,#MAPTYPE_TERRAIN)
   SaveImage(planetSide(i)\numTexture,"temp"+Str(i)+".bmp")
   FreeImage(planetSide(i)\numTexture)
  
  planetSide(i)\numTexture = LoadTexture(#PB_Any,"temp"+Str(i)+".bmp")
  DeleteFile("temp"+Str(i)+".bmp")
  
  planetSide(i)\numMaterial = CreateMaterial(#PB_Any,TextureID(planetSide(i)\numTexture))            
  MaterialFilteringMode(planetSide(i)\numMaterial, #PB_Material_None)
Next i

;- Mesh and entity
planetEntity = createPlanet(planetSide(),tesselation)

; Attach a node, so you can add other elements later (atmosphere, moons, etc..)
CreateNode(0,0,0,0)
AttachNodeObject(0,EntityID(planetEntity),#PB_Node_Entity)
  

;- Camera 
CreateCamera(0, 0, 0 , 100 , 100) 
MoveCamera(0,0,1000,1500) 
CameraLookAt(0,0,0,0) 

;- Light 
AmbientColor($333333) 
CreateLight(0,$FFFFFF,250,500,1000)


;- Main loop
Repeat 
  While WindowEvent():Wend
  Delay(1)
   
  ;- Mouse
  ExamineMouse()
  anglex+MouseDeltaX()
  angley+MouseDeltaY()
  movez = MouseWheel()

  ;- Keyboard
  ExamineKeyboard() 
  
  If KeyboardReleased(#PB_Key_W)
    wireMode = 1-wireMode
    If wireMode = 1
      CameraRenderMode(0,#PB_Camera_Wireframe)
    Else
      CameraRenderMode(0,#PB_Camera_Textured)
    EndIf 
  EndIf
  
  If KeyboardReleased(#PB_Key_R) And tesselation > 1
    DetachNodeObject(0,EntityID(planetEntity),#PB_Node_Entity)
    FreeEntity(planetEntity)
    tesselation-1
    Debug tesselation
    planetEntity = createPlanet(planetSide(),tesselation)
    AttachNodeObject(0,EntityID(planetEntity),#PB_Node_Entity)
  EndIf
  
  If KeyboardReleased(#PB_Key_T) And tesselation < 10
    DetachNodeObject(0,EntityID(planetEntity),#PB_Node_Entity)
    FreeEntity(planetEntity)
    tesselation+1
    planetEntity = createPlanet(planetSide(),tesselation)
    AttachNodeObject(0,EntityID(planetEntity),#PB_Node_Entity)
  EndIf
  
  If KeyboardPushed(#PB_Key_Right)
    RotateNode(0,0,-1,0,#PB_Relative)
    autoRotate = #False
  EndIf
  If KeyboardPushed(#PB_Key_Left)
    RotateNode(0,0,1,0,#PB_Relative)
    autoRotate = #False
  EndIf
 
  If KeyboardPushed(#PB_Key_Up)
    RotateNode(0,1,0,0,#PB_Relative)
    autoRotate = #False
  EndIf
  If KeyboardPushed(#PB_Key_Down)
    RotateNode(0,-1,0,0,#PB_Relative)
    autoRotate = #False
  EndIf
 
  
  MoveCamera(0,0,0,movez * -50)
  
  ; show it all
  RenderWorld() 
  
  ; Flip buffers to avoid tearing  
  FlipBuffers()
Until KeyboardPushed(#PB_Key_Escape) 
End


;************************************************************************************
;-                                 ---- Datas ----
;************************************************************************************
DataSection
  terrainTypes:
  Data.s "Tundra"
  Data.f 1,4,0,2.5,0,1.5
  Data.i $79A700
  
  Data.s "Grass"
  Data.f 1,3,0,2,18,3.5
  Data.i $5AD00F
  
  Data.s "Forest" 
  Data.f 1.1,3,5,3,25,2
  Data.i $3EA300
  
  Data.s "Desert" 
  Data.f 1.3,5,0,1.25,50,2
  Data.i $E5DFAD
  
  Data.s "Mountain"
  Data.f 1.7,2,20,3,10,2
  Data.i $BBBBBB
  
  Data.s "XXX END XXX"
  
  meshDatas:
  ; Nb vertices, nb triangles
  Data.i 4,2
  ; Vertices: pos / normals / uv
  Data.f -0.5,0.5,-0.5
  Data.f 0,1,0
  Data.f 0,0
  Data.f 0.5,0.5,-0.5
  Data.f 0,1,0
  Data.f 1,0
  Data.f 0.5,0.5,0.5
  Data.f 0,1,0
  Data.f 1,1
  Data.f -0.5,0.5,0.5
  Data.f 0,1,0
  Data.f 0,1
  ; Faces
  Data.i 2,1,0
  Data.i 0,3,2

EndDataSection
applePi
Addict
Addict
Posts: 1404
Joined: Sun Jun 25, 2006 7:28 pm

Re: [4.60] Telluric Planet Generator

Post by applePi »

this is a great demo, tasty planet, its texture like a velvet clothing
, or like an alien cartoon planet.
and a code resource for a year. i was looking at the Venus and Jupiter conjunction - interactive when i read your post here about this planet.
thanks
User avatar
J@ckWhiteIII
Enthusiast
Enthusiast
Posts: 183
Joined: Fri May 25, 2012 7:39 pm

Re: [4.60] Telluric Planet Generator

Post by J@ckWhiteIII »

Hm..launching it works fine for me, except for one thing: all I get is a black screen in the window. the debug stuff seems to be fine
Post Reply