FreeBASIC マニュアルのトップに戻る

FreeBASIC terminal velocity sim

目次→フォーラム→FreeBASIC→補足Ray Tracer←オリジナル・サイト

光投射(光と影) 左にメニュー・フレームが表示されていない場合は、ここをクリックして下さい

←リンク元に戻る プログラム開発関連に戻る

光投射(光と影)  UEZ さんが JavaScript コードから FB に変換したものに、D.J.Peters さんがベクトル数学をマクロに置き換えて処理を高速化したものです。

機能の一覧:
1)独自の位置、色、サイズを持つ複数の球
2)球体の最終的な色は拡散、反射、鏡面反射の色から組み合わされています
3)織柄付きの床 - タイル張りのカーペット
4)床は世界の光沢のある反射を持っています
5)空はグラデーションと「星」でテクスチャリングされています
6)床と球の柔らかい影
7)被写界深度も存在します

'Ray Tracer
'by D.J.Peters ≫ Sep 21, 2022 22:14
'https://www.freebasic.net/forum/viewtopic.php?p=294787#p294787


#Ifndef __FB_64BIT__
 #cmdline "-gen gcc -arch pentium4-sse3 -Wc -O3 -fpu sse -O 3 -fpmode fast -asm intel"
 Type real As Single
#Else
 #cmdline "-arch x86-64 -Wc -O3 -fpmode fast -fpu sse -O 3 -asm intel"
 Type real As Double
#EndIf

' Ported from https://js1k.com/2017-magic/demo/2648 by Igor Sbitnev To FB by UEZ build 2021-01-19

' original on old laptop ~24 seconds optimized by d.j.peters ~5 seconds now :-)

#Include "crt.bi"

Type Vector3D
  As real x, y, z
End Type

Type tSpheresData
   As Vector3D position, colour
   As real radius,r2
End Type

Enum eIntersection 
  INTERSECTION_NONE
  INTERSECTION_SPHERE
  INTERSECTION_FLOOR
End Enum  

Type tIntersect
   As eIntersection intersection
   As Vector3D rayEnd, lightDir, normal, sphereColor
End Type

#define vector(_x,_y,_z) Type<vector3d>(_x, _y, _z)
#define vCross(_a,_b) vector(_a.y*_b.z - _a.z*_b.y, _a.z*_b.x - _a.x*_b.z, _a.x*_b.y - _a.y*_b.x)
Function vNormal(ByRef v As Const Vector3D) As Vector3D
  Dim As real l = v.x*v.x + v.y*v.y + v.z*v.z
  If l Then l=1/Sqr(l)
  Return vector(v.x*l,v.y*l,v.z*l)
End Function

Dim Shared As Vector3D CENTER_TILE_COLOR
Dim Shared As Vector3D OTHERS_TILE_COLOR
Dim Shared As Vector3D COLOR_SKY
Dim Shared As Vector3D COLOR_LIGHT_SOURCE
Dim Shared As Vector3D COLOR_STARS
CENTER_TILE_COLOR  = Vector(8, 0, 8)
OTHERS_TILE_COLOR  = Vector(8, 5, 8)
COLOR_SKY          = Vector(5, 6, 8)
COLOR_LIGHT_SOURCE = Vector(8, 8, 8)
COLOR_STARS        = Vector(4, 4, 8)

Dim Shared As tSpheresData spheresData(5)
                                 ' position ,color,radius, radius squared
spheresData(0)=Type<tSpheresData>(Vector(10.0,  2.0, 2.0),Vector(4, 0, 4),1.5,1.5*1.5)
spheresData(1)=Type<tSpheresData>(Vector(-3.0,  0.0, 2.0),Vector(8, 5, 7),1.5,1.5*1.5)
spheresData(2)=Type<tSpheresData>(Vector( 3.0,  0.0, 2.0),Vector(0, 0, 4),1.5,1.5*1.5)
spheresData(3)=Type<tSpheresData>(Vector( 1.5,  0.0, 4.5),Vector(8, 8, 6),1.5,1.5*1.5)
spheresData(4)=Type<tSpheresData>(Vector(-1.0, 10.0, 4.0),Vector(0, 4, 4),4.0,4.0*4.0)
spheresData(5)=Type<tSpheresData>(Vector( 0.0,  0.0, 7.0),Vector(8, 5, 4),1.5,1.5*1.5)



Const RAYS_PER_PIXEL = 32
Dim Shared As real preX(RAYS_PER_PIXEL-1)
Dim Shared As real preY(RAYS_PER_PIXEL-1)

Function trace(ByRef rayStart As Const Vector3D, _
               ByRef rayDir   As Const Vector3D) As tIntersect
   Dim As Vector3D rayEnd=Any
   Dim As Vector3D lightDir=Any
   Dim As Vector3D normal=Any
   Dim As Vector3D sphereColor=Any
   Dim As Vector3D raySphere=Any
   Dim As Vector3D tmp3=Any
   Dim As Vector3D lightPos = Vector(Rnd() * 27, -81 + Rnd() * 27, 81)
   Dim As real b=Any, c=Any, d=Any, tSphere=Any,tmp=Any
   Dim As real tFloor = -rayStart.z / rayDir.z
   Dim As real tNear = 2^31
   Dim As eIntersection intersection = INTERSECTION_NONE
   Dim As Integer sphereIndex=-1
   If tFloor > 0 And rayStart.z > 0 Then
      tNear=tFloor
   End If
   For i As Integer = 0 To UBound(spheresData)
      'raySphere = vSub(spheresData(i).position,rayStart)
      With spheresData(i).position  
        raySphere.x = .x-rayStart.x
        raySphere.y = .y-rayStart.y
        raySphere.z = .z-rayStart.z
      End With  
      b = raySphere.x*rayDir.x _
        + raySphere.y*rayDir.y _
        + raySphere.z*rayDir.z
      If b < 0 Then Continue For
      c = raySphere.x*raySphere.x _
        + raySphere.y*raySphere.y _
        + raySphere.z*raySphere.z
      c-= b*b  
      d = spheresData(i).r2 
      If c > d Then Continue For
      tSphere = b - Sqr(d - c)
      If tSphere<0 OrElse tSphere > tNear Then Continue For
      sphereIndex = i : tNear = tSphere
    
   Next
   If sphereIndex>-1 Then
     intersection = INTERSECTION_SPHERE
     'sphereColor = spheresData(sphereIndex).colour
     memcpy(@sphereColor.x,@spheresData(sphereIndex).colour.x,SizeOf(real)*3)
     rayEnd.x = rayStart.x + rayDir.x*tNear
     rayEnd.y = rayStart.y + rayDir.y*tNear
     rayEnd.z = rayStart.z + rayDir.z*tNear
     'lightDir = vNormal(vSub(lightPos,rayEnd))
     lightDir.x = lightPos.x - rayEnd.x
     lightDir.y = lightPos.y - rayEnd.y
     lightDir.z = lightPos.z - rayEnd.z
     tmp = lightDir.x*lightDir.x
     tmp+= lightDir.y*lightDir.y
     tmp+= lightDir.z*lightDir.z
     If tmp Then
       tmp=1/Sqr(tmp)
       lightDir.x*=tmp
       lightDir.y*=tmp
       lightDir.z*=tmp
     EndIf  
     'normal = vNormal(vSub(rayEnd,spheresData(sphereIndex).position))
     normal.x = rayEnd.x - spheresData(sphereIndex).position.x
     normal.y = rayEnd.y - spheresData(sphereIndex).position.y
     normal.z = rayEnd.z - spheresData(sphereIndex).position.z
     tmp = normal.x*normal.x
     tmp+= normal.y*normal.y
     tmp+= normal.z*normal.z
     If tmp Then
       tmp=1/Sqr(tmp)
       normal.x*=tmp
       normal.y*=tmp
       normal.z*=tmp
     EndIf  
     
   ElseIf tFloor > 0 And rayStart.z > 0 Then
     intersection = INTERSECTION_FLOOR
     'rayEnd = vAdd(rayStart,vScale(rayDir,tFloor))
     rayEnd.x = rayStart.x + rayDir.x*tFloor
     rayEnd.y = rayStart.y + rayDir.y*tFloor
     rayEnd.z = rayStart.z + rayDir.z*tFloor
     'lightDir = vNormal(vSub(lightPos, rayEnd))
     lightDir.x = lightPos.x - rayEnd.x
     lightDir.y = lightPos.y - rayEnd.y
     lightDir.z = lightPos.z - rayEnd.z
     tmp = lightDir.x*lightDir.x
     tmp+= lightDir.y*lightDir.y
     tmp+= lightDir.z*lightDir.z
     If tmp Then
       tmp=1/Sqr(tmp)
       lightDir.x*=tmp
       lightDir.y*=tmp
       lightDir.z*=tmp
     EndIf  
     normal = Vector(0, 0, 1)
   EndIf
   Return Type<tIntersect>(intersection, rayEnd, lightDir, normal, sphereColor)
End Function

Function getFloorColor(x As Integer, y As Integer) As Vector3D
   If x + y <> 0 Then
      If fmod(x, 3) = 1 And fmod(y, 3) = 1 Then Return CENTER_TILE_COLOR
      Return getFloorColor((x \ 3), (y \ 3))
   End If
   Return OTHERS_TILE_COLOR
End Function

Function sample(rayStart As Vector3D, rayDir As Vector3D, renderStars As Boolean = FALSE) As Vector3D
  Static As Integer j=0
  Dim As Vector3D refDir=Any
  Dim As Vector3D colour=Any
  Dim As Vector3D diffuse=Any
  Dim As Vector3D specular=Any
  Dim As Vector3D reflection=Any
  Dim As Vector3D floorColor=Any
  Dim As Vector3D uVector=Any
  Dim As Vector3D vVector=Any
  Dim As Vector3D rndDir=Any
  Dim As Vector3D tmp3=Any
  Dim As real tmp = Any,tmp2 
  Dim As tIntersect ti = trace(rayStart, rayDir)
  Select Case As Const ti.intersection
  Case INTERSECTION_NONE
    If renderStars AndAlso Rnd() > 0.5 Then Return COLOR_STARS
    tmp = Pow(1 - rayDir.z, 4)
    tmp3.x = COLOR_SKY.x*tmp
    tmp3.y = COLOR_SKY.y*tmp
    tmp3.z = COLOR_SKY.z*tmp
    Return tmp3'Iif(renderStars And Rnd() > 0.5, COLOR_STARS, vScale(COLOR_SKY,tmp))
  Case INTERSECTION_SPHERE
    'tmp = vDot(ti.normal, ti.lightDir)*0.7
    tmp = ti.normal.x*ti.lightDir.x
    tmp+= ti.normal.y*ti.lightDir.y
    tmp+= ti.normal.z*ti.lightDir.z
    tmp*=0.7
    'diffuse = vScale(ti.sphereColor,tmp)
    colour.x = ti.sphereColor.x*tmp
    colour.y = ti.sphereColor.y*tmp
    colour.z = ti.sphereColor.z*tmp
    'tmp=Pow(vDot(ti.normal, vNormal(vSub(ti.lightDir, rayDir))), 64)
    tmp3.x = ti.lightDir.x - rayDir.x
    tmp3.y = ti.lightDir.y - rayDir.y
    tmp3.z = ti.lightDir.z - rayDir.z
    tmp = tmp3.x*tmp3.x
    tmp+= tmp3.y*tmp3.y
    tmp+= tmp3.z*tmp3.z
    If tmp Then
      tmp=1/Sqr(tmp)
      tmp3.x*=tmp
      tmp3.y*=tmp
      tmp3.z*=tmp
    EndIf
    tmp = ti.normal.x*tmp3.x
    tmp+= ti.normal.y*tmp3.y
    tmp+= ti.normal.z*tmp3.z
    tmp = pow(tmp,64)
    ' specular = vScale(COLOR_LIGHT_SOURCE,tmp)
    colour.x += COLOR_LIGHT_SOURCE.x*tmp
    colour.y += COLOR_LIGHT_SOURCE.y*tmp
    colour.z += COLOR_LIGHT_SOURCE.z*tmp
    ' tmp = -2 * vDot(ti.normal, rayDir)
    tmp = ti.normal.x*rayDir.x
    tmp+= ti.normal.y*rayDir.y
    tmp+= ti.normal.z*rayDir.z
    tmp*=-2
    ' refDir = vAdd(rayDir,vScale(ti.normal,tmp))
    refDir.x = rayDir.x + ti.normal.x*tmp
    refDir.y = rayDir.y + ti.normal.y*tmp
    refDir.z = rayDir.z + ti.normal.z*tmp
    tmp3=sample(ti.rayEnd, refDir)
    'reflection = vScale(tmp3,0.4)
    reflection.x = tmp3.x*0.4
    reflection.y = tmp3.y*0.4
    reflection.z = tmp3.z*0.4
    'colour = vAdd(diffuse,vAdd(specular,reflection))
    colour.x +=  reflection.x
    colour.y +=  reflection.y
    colour.z +=  reflection.z
  Case INTERSECTION_FLOOR
    'tmp = -2 * vDot(ti.normal, rayDir)
    tmp = ti.normal.x*rayDir.x
    tmp+= ti.normal.y*rayDir.y
    tmp+= ti.normal.z*rayDir.z
    tmp*=-2
    'refDir = vAdd(rayDir,vScale(ti.normal,tmp))
    refDir.x = rayDir.x + ti.normal.x*tmp
    refDir.y = rayDir.y + ti.normal.y*tmp
    refDir.z = rayDir.z + ti.normal.z*tmp
    floorColor = getFloorColor(CInt(fmod((ti.rayEnd.x + 81) * 27, 81)), CInt(fmod((ti.rayEnd.y + 81) * 27, 81)))
    ' uVector = vCross(rayDir,refDir)
    uVector.x = rayDir.y*refDir.z - rayDir.z*refDir.y
    uVector.y = rayDir.z*refDir.x - rayDir.x*refDir.z
    uVector.z = rayDir.x*refDir.y - rayDir.y*refDir.x
    ' vVector = vCross(uVector,refDir)
    vVector.x = uVector.y*refDir.z - uVector.z*refDir.y
    vVector.y = uVector.z*refDir.x - uVector.x*refDir.z
    vVector.z = uVector.x*refDir.y - uVector.y*refDir.x
    'rndDir = vAdd(refDir,vAdd(vScale(uVector,preX(j)),vScale(vVector,preY(j))))    
    tmp = preX(j) : tmp2 = preY(j)
    rndDir.x = refDir.x + uVector.x*tmp + vVector.x*tmp2
    rndDir.y = refDir.y + uVector.y*tmp + vVector.y*tmp2
    rndDir.z = refDir.z + uVector.z*tmp + vVector.z*tmp2
    tmp3 = sample(ti.rayEnd,rndDir)
    'colour = vScale(vAdd(floorColor,tmp3),.5)
    colour.x = floorColor.x + tmp3.x : colour.x*=.5
    colour.y = floorColor.y + tmp3.y : colour.y*=.5
    colour.z = floorColor.z + tmp3.z : colour.z*=.5
    j+=1 : If j=RAYS_PER_PIXEL Then j=0
  End Select
  If trace(ti.rayEnd, ti.lightDir).intersection Then 
    'colour = vScale(colour,.5)
    colour.x *= .5
    colour.y *= .5
    colour.z *= .5
  EndIf  
  Return colour
End Function



Const As Integer iW = 512
Const As Integer iH = 512
Const DISTANCE_TO_VIEWPORT = 10
Const VIEWPORT_WIDTH = 12
Const VIEWPORT_HEIGHT = 12
Const ALPHA_CHANNEL_COLOR = 255
Const As real scaleX = VIEWPORT_WIDTH / iW
Const As real ScaleY = VIEWPORT_HEIGHT / iH


Dim As Vector3D rayStart=Any
Dim As Vector3D rayDir=Any
Dim As Vector3D viewportPixel=Any
Dim As Vector3D viewportCenter = Any
Dim As Vector3D leftDown = Any
Dim As Vector3D UP_DIRECTION = Vector(0, 0, 1)
Dim As Vector3D camera = Vector(-7, -10, 8)
Dim As Vector3D target = Vector(0, 0, 4)
Dim As Vector3D normalToViewport = vNormal(vector(camera.x-target.x,camera.y-target.y,camera.z-target.z ))
Dim As Vector3D uVector = vNormal(vCross(UP_DIRECTION, normalToViewport))
Dim As Vector3D vVector = vCross(uVector, normalToViewport)

viewportCenter.x = camera.x + normalToViewport.x*-DISTANCE_TO_VIEWPORT
viewportCenter.y = camera.y + normalToViewport.y*-DISTANCE_TO_VIEWPORT
viewportCenter.z = camera.z + normalToViewport.z*-DISTANCE_TO_VIEWPORT

leftDown.x = viewportCenter.x + uVector.x*-VIEWPORT_WIDTH/2 + vVector.x*-VIEWPORT_HEIGHT/2
leftDown.y = viewportCenter.y + uVector.y*-VIEWPORT_WIDTH/2 + vVector.y*-VIEWPORT_HEIGHT/2
leftDown.z = viewportCenter.z + uVector.z*-VIEWPORT_WIDTH/2 + vVector.z*-VIEWPORT_HEIGHT/2

ScreenRes iW, iH, 32,2
ScreenSet 1,0

#Macro colorClamp(_c,_v) 
  If _v>255 Then
   _c=255
  ElseIf _v<0 Then
   _c=0
  Else
   _c=_v
  EndIf  
#EndMacro  

Dim As Any Ptr img = ImageCreate(iW,1)
Dim As ULong Ptr pixels,row
Dim As Long pitch
ImageInfo img,,,,pitch,pixels 
pitch Shr=2 ' <-- from bytes to ulong
Dim As ULong r,g,b

Dim As Vector3d preScaleX(iW)
For x As Integer=0 To iW-1
  Dim As real tmp=x*ScaleX
  preScaleX(x).x = uVector.x*tmp
  preScaleX(x).y = uVector.y*tmp
  preScaleX(x).z = uVector.z*tmp
Next
For i As Integer=0 To RAYS_PER_PIXEL-1
  preX(i)=(Rnd() - 0.5) / 3
  preY(i)=(Rnd() - 0.5) / 3
Next
Dim As real tmp,tmp2,onePercent = 100/iH
Dim As Integer oldPercent,newPercent
Dim As Vector3d tmp3,preScaleY
Dim As Double tAll = Timer()
For y As Integer = 0 To iH-1
  row = pixels 
  tmp=y*ScaleY
  preScaleY.x = vVector.x*tmp
  preScaleY.y = vVector.y*tmp
  preScaleY.z = vVector.z*tmp  
  For x As Integer = 0 To iW-1
    Dim As Vector3D colorSum = Vector(0, 0, 0)
    For j As Integer = 0 To RAYS_PER_PIXEL - 1
      'rayStart = vAdd(camera,vAdd(vScale(uVector,preX(j)),vScale(vVector,preY(j))))
      tmp = preX(j) : tmp2 = preY(j)
      rayStart.x = camera.x + uVector.x*tmp + vVector.x*tmp2
      rayStart.y = camera.y + uVector.y*tmp + vVector.y*tmp2
      rayStart.z = camera.z + uVector.z*tmp + vVector.z*tmp2
      'viewportPixel = vSub(vAdd(leftDown,vAdd(preScaleX(x),preScaleY)),rayStart)
      tmp3 = preScaleX(x)
      rayDir.x = leftDown.x + tmp3.x + preScaleY.x : rayDir.x-=rayStart.x
      rayDir.y = leftDown.y + tmp3.y + preScaleY.y : rayDir.y-=rayStart.y
      rayDir.z = leftDown.z + tmp3.z + preScaleY.z : rayDir.z-=rayStart.z
      tmp = rayDir.x*rayDir.x
      tmp+= rayDir.y*rayDir.y
      tmp+= rayDir.z*rayDir.z
      If tmp Then
        tmp=1/Sqr(tmp)
        rayDir.x *= tmp
        rayDir.y *= tmp
        rayDir.z *= tmp
      EndIf
      'rayDir = vNormal(viewportPixel)
      'rayDir.x = viewportPixel.x
      'rayDir.y = viewportPixel.y
      'rayDir.z = viewportPixel.z
      tmp3=sample(rayStart, rayDir, TRUE)
      colorSum.x+= tmp3.x
      colorSum.y+= tmp3.y
      colorSum.z+= tmp3.z
    Next
    colorClamp(r,colorSum.x)
    colorClamp(g,colorSum.y)
    colorClamp(b,colorSum.z)
    row[x]=RGB(r,g,b)
  Next
  row+=pitch
  Put (0,y),img,PSet
  Flip 'screenunlock
  newPercent = y*onePercent
  If newPercent<>oldPercent Then
    WindowTitle "done: " & newPercent & " %"
    oldPercent=newPercent
  EndIf  
  If ((y Mod 8)=0) AndAlso Inkey()<>"" Then Exit For
Next
tAll = Timer()-tAll
WindowTitle("Ray Tracer / Rendered in " & Int(tAll) & " seconds done ...")
Sleep
 
補足 に戻る
←リンク元に戻る プログラム開発関連に戻る
ページ歴史:2022-09-29
日本語翻訳:WATANABE Makoto、原文著作者:UEZ、D.J.Peters

ホームページのトップに戻る

表示-非営利-継承