ピクセルシェーダで空を描く
![f:id:agehama:20150623213441g:plain f:id:agehama:20150623213441g:plain](http://cdn-ak.f.st-hatena.com/images/fotolife/a/agehama/20150623/20150623213441.gif)
Siv3Dが(2D限定で)ピクセルシェーダをサポートしたらしいので、GPU Gems 2に載っていた大気散乱シミュレーションのコードを移植してみました。頂点シェーダからそのままピクセルシェーダに移しただけなので少し遅いです。
大気散乱シミュレーションについては、このページの参考文献に高品質な実装法の紹介があるので、興味のある人はそちらを読むといいと思います。(自分はインクを塗るのに忙しくてまだ読んでない。)
あと、地形生成のソースコードはここからお借りしました。
以下スクリーンショットとソースコードです。
ソースコード
PS_Sky.hlsl
cbuffer CB : register(b1) { float3 frustumRayTL; float3 frustumRayTR; float3 frustumRayBL; uint2 windowSize; }; cbuffer SkyCB : register(b2) { float3 v3CameraPos; // The camera's current position float fCameraHeight; // The camera's current height float3 v3LightPos; // The direction vector to the light source float fCameraHeight2; // fCameraHeight^2 float3 v3InvWavelength; // 1 / pow(wavelength, 4) for the red, green, and blue channels float fScale; // 1 / (fOuterRadius - fInnerRadius) float fOuterRadius; // The outer (atmosphere) radius float fOuterRadius2; // fOuterRadius^2 float fInnerRadius; // The inner (planetary) radius float fInnerRadius2; // fInnerRadius^2 float fKrESun; // Kr * ESun float fKmESun; // Km * ESun float fKr4PI; // Kr * 4 * PI float fKm4PI; // Km * 4 * PI float fScaleDepth; // The scale depth (i.e. the altitude at which the atmosphere's average density is found) float fScaleOverScaleDepth; // fScale / fScaleDepth float g; float exposure; }; struct VS_OUTPUT { float4 postion : SV_POSITION; }; float3 IntersectionPos(float3 rayPos, float3 rayDir, float sphereRadius) { const float A = dot(rayDir, rayDir); const float B = 2.0*dot(rayPos, rayDir); const float C = dot(rayPos, rayPos) - sphereRadius*sphereRadius; return B*B - 4.0*A*C < 0.000001 ? float3(0, 0, 0) : (rayPos + rayDir*(0.5*(-B + sqrt(B*B - 4.0*A*C)) / A)); } float IntegralApproximation(float fCos) { float x = 1.0 - fCos; return fScaleDepth * exp(-0.00287 + x*(0.459 + x*(3.83 + x*(-6.80 + x*5.25)))); } float4 PS(in VS_OUTPUT input) : SV_Target { float2 uv = (input.postion.xy / float2(windowSize)); float3 deltaX = (frustumRayTR - frustumRayTL)*uv.x; float3 deltaY = (frustumRayBL - frustumRayTL)*uv.y; float3 ray = normalize((frustumRayTL + deltaX + deltaY).xyz); float3 skyPos = IntersectionPos(v3CameraPos.xyz, ray, fOuterRadius); if (skyPos.x == 0) { //Camera is out of celestial sphere. clip(-1); } float3 invWavelength = v3InvWavelength; float3 v3Pos = skyPos; float3 v3Ray = v3Pos - v3CameraPos.xyz; float fFar = length(v3Ray); v3Ray /= fFar; // Calculate the ray's starting position, then calculate its scattering offset float3 v3Start = v3CameraPos; float fHeight = length(v3Start); float fDepth = exp(fScaleOverScaleDepth * (fInnerRadius - fCameraHeight)); float fStartAngle = dot(v3Ray, v3Start) / fHeight; float fStartOffset = fDepth*IntegralApproximation(fStartAngle); const int nSamples = 5; const float fSamples = nSamples; // Initialize the scattering loop variables float fSampleLength = fFar / fSamples; float fScaledLength = fSampleLength * fScale; float3 v3SampleRay = v3Ray * fSampleLength; float3 v3SamplePoint = v3Start + v3SampleRay * 0.5; // Now loop through the sample rays float3 v3FrontColor = float3(0.0, 0.0, 0.0); for (int i = 0; i<nSamples; ++i) { float fHeight = length(v3SamplePoint); float fDepth = exp(fScaleOverScaleDepth * (fInnerRadius - fHeight)); float fLightAngle = dot(v3LightPos, v3SamplePoint) / fHeight; float fCameraAngle = dot(v3Ray, v3SamplePoint) / fHeight; float fScatter = (fStartOffset + fDepth*(IntegralApproximation(fLightAngle) - IntegralApproximation(fCameraAngle))); float3 v3Attenuate = exp(-fScatter * (invWavelength * fKr4PI + fKm4PI)); v3FrontColor += v3Attenuate * (fDepth * fScaledLength); v3SamplePoint += v3SampleRay; } float4 primaryColor = float4(v3FrontColor * (invWavelength * fKrESun), 1.0); float4 secondaryColor = float4(v3FrontColor * fKmESun, 1.0); float3 v3Direction = v3CameraPos - v3Pos; float fCos = dot(v3LightPos, v3Direction) / length(v3Direction); float fRayPhase = 0.75 * (1.0 + fCos * fCos); const float g2 = g*g; float fMiePhase = 1.5 * ((1.0 - g2) / (2.0 + g2)) * (1.0 + fCos*fCos) / pow(abs(1.0 + g2 - 2.0*g*fCos), 1.5); float3 raycolor = (primaryColor*fRayPhase).xyz; float3 miecolor = (secondaryColor*fMiePhase).xyz; float3 c = float3(1.0, 1.0, 1.0) - exp(-exposure * (raycolor + miecolor)); return float4(saturate(c), 1.0); }
Main.cpp
#include <Siv3D.hpp> //方向を操作するためのGUI class GuideArrow { public: GuideArrow() : transform(Mat4x4::Identity()), cone(MeshData::Cone(0.25)), cylinder(MeshData::Cylinder(0.1, 1.5)), sphere(MeshData::Sphere(0.2)) {} //矢印の描画位置は常にカメラから見て同じ位置にあるべきなので、カメラの位置情報からtransformを更新する void update(double rotateSpeed = 5.0) { const Vec3 newPosition = Graphics3D::GetCamera().pos + Graphics3D::ToRay(Window::Center()).direction*5.0; Quaternion currentRotation = rotation(); transform = Mat4x4::Rotate(currentRotation).translated(newPosition); if (Input::MouseR.pressed && !Mouse::Delta().isZero) { Vec2 delta = Mouse::Delta(); const Vec3 previousMouseRay = Graphics3D::ToRay(Mouse::PreviousPos()).direction.normalized(); const Vec3 currentMouseRay = Graphics3D::ToRay(Mouse::Pos()).direction.normalized(); const Vec3 normal = previousMouseRay.cross(currentMouseRay).normalized(); const auto deltaRotation = s3d::Quaternion(normal, Radians(-rotateSpeed)); //マウスドラッグに合わせて矢印を回転する transform = Mat4x4::Rotate(currentRotation).rotated(deltaRotation).translated(newPosition); } } void drawForward(double scale = 1.0, const Vec3& move = { 0, 0, 0 }) { const Mat4x4 drawTransform = Mat4x4::Scale(scale)*transform.translated(move); cone.drawForward(Mat4x4::Translate(Float3{ 0, 1.35, 0 })*drawTransform); cylinder.drawForward(Mat4x4::Translate(Float3{ 0, 0.75, 0 })*drawTransform); sphere.drawForward(Mat4x4::Translate(Float3{ 0, 0, 0 })*drawTransform); } //矢印の向いている方向を取得する Vec3 getDirection()const { const Vector originalArrowDirection = s3d::ToVector(Graphics3D::GetCamera().up, 0.0); const Vector rotatedArrowDirection = XMVector3Transform(originalArrowDirection, Mat4x4::Rotate(rotation())); const Vector& v = rotatedArrowDirection; return Vec3(DirectX::XMVectorGetX(v), DirectX::XMVectorGetY(v), DirectX::XMVectorGetZ(v)); } void setDirection(const Vec3& direction) { const Vec3 newPosition = Graphics3D::GetCamera().pos + Graphics3D::ToRay(Window::Center()).direction*5.0; const Vec3 original = Vec3(0, 1, 0); const double angle = Acos(original.dot(direction)); transform = Mat4x4::Rotate(Quaternion(original.cross(direction).normalized(), angle)).translated(newPosition); } void setRotation(const Quaternion& quaternion) { const Vec3 newPosition = Graphics3D::GetCamera().pos + Graphics3D::ToRay(Window::Center()).direction*5.0; transform = Mat4x4::Rotate(quaternion).translated(newPosition); } void reset() { setDirection(Vec3(0, 1, 0)); } private: Quaternion rotation()const { Vec3 position; Vec3 scale; Quaternion rotation; if (!transform.decompose(scale, rotation, position)) { LOG_ERROR(L"行列の分解に失敗しました。"); return Quaternion::Identity(); } return rotation; } Mat4x4 transform; Mesh cone; Mesh cylinder; Mesh sphere; }; Vec3 SunDirection(int day, double time, double eastLongitude, double northLatitude); class PS_Sky { public: PS_Sky(): shader(L"PS_Sky.hlsl"), gui(GUIStyle::Default), guiActiveness(true) { gui.setTitle(L"大気散乱のシミュレーション"); const int widthSlider = 256; const float widthConstant = 0.005f; gui.add(GUIText::Create(L"レイリー散乱定数:\t")); gui.addln(L"Kr", GUISlider::Create(0.0025f - widthConstant, 0.0025f + widthConstant, 0.0025f, widthSlider)); gui.add(GUIText::Create(L"ミー散乱定数:\t\t")); gui.addln(L"Km", GUISlider::Create(0.0010f - widthConstant, 0.0010f + widthConstant, 0.0010f, widthSlider)); gui.add(GUIText::Create(L"明るさ:\t\t\t")); gui.addln(L"Exposure", GUISlider::Create(0.001, 0.3, 0.05, widthSlider)); gui.add(GUIText::Create(L"光の方向操作:\t\t")); gui.addln(L"Auto", GUIRadioButton::Create({ L"太陽", L"サイン波", L"手動([Alt]キーで減速)" }, 0u)); gui.add(L"Save", GUIButton::Create(L"画像を保存")); gui.addln(L"Reset", GUIButton::Create(L"全てリセット")); gui.add(GUIHorizontalLine::Create()); gui.add(GUIText::Create(L"WASDキーと矢印キーでカメラの移動と回転。\n右クリック&ドラッグで光の方向を操作。")); gui.setPos(1280, 370); directionTimer.start(); } void update() { if (Input::KeyH.clicked) { guiActiveness = !guiActiveness; guiActiveness ? gui.show() : gui.hide(); } //スクリーンショットを撮り終わったら再びUIを表示する if (ScreenCapture::HasNewFrame()) { guiActiveness ? gui.show() : gui.hide(); } updateLightDirection(); updateCB(); updateCBSky(directionGUI.getDirection()); Graphics2D::BeginShader(shader); Graphics2D::SetConstant(ShaderStage::Pixel, 1, cb); Graphics2D::SetConstant(ShaderStage::Pixel, 2, cbSky); Window::ClientRect().draw(); Graphics2D::EndShader(); Graphics::Render2DBackground(); //スクリーンショットを撮る時は一旦UIを隠す if (gui.button(L"Save").pushed) { gui.hide(); ScreenCapture::Save(); } else if (guiActiveness) { directionGUI.drawForward(0.5, -1.0*Graphics3D::GetCamera().up); } } void reset() { guiActiveness = true; gui.show(); gui.setPos(1280, 370); gui.slider(L"Kr").setValue(0.0025f); gui.slider(L"Km").setValue(0.0010f); gui.slider(L"Exposure").setValue(0.05); gui.radioButton(L"Auto").check(0u); directionGUI.reset(); directionTimer.restart(); } Vec3 sunDirection()const { return directionGUI.getDirection(); } bool pushedResetButton()const { return gui.button(L"Reset").pushed; } private: void updateLightDirection() { if (guiActiveness && Input::MouseR.pressed) { gui.radioButton(L"Auto").check(2u); } if (gui.radioButton(L"Auto").hasChanged) { directionTimer.restart(); } if (gui.radioButton(L"Auto").checkedItem) { const unsigned key = gui.radioButton(L"Auto").checkedItem.value(); //太陽 if (key == 0u) { const int T = 10000; const double time = 24.0*(directionTimer.elapsed() % T) / T; const int day = (directionTimer.elapsed() % (T * 365)) / T; directionGUI.setDirection(SunDirection(day, time, 139.77, 35.68)); } //サイン波 else if (key == 1u) { const int T_XZ = 10000; const double rate_XZ = 1.0*(directionTimer.elapsed() % T_XZ) / T_XZ; const int T_Y = 1000; const double rate_Y = 1.0*(directionTimer.elapsed() % T_Y) / T_Y; const Vec3 direction(Cos(rate_XZ*TwoPi), 0.25*Sin(rate_Y*TwoPi), Sin(rate_XZ*TwoPi)); directionGUI.setDirection(direction.normalized()); } //手動 else { directionGUI.update(Input::KeyAlt.pressed ? 0.5 : 5.0); } } } void updateCB() { const auto window = Window::ClientRect(); cb->frustumRayTL = Float4(Graphics3D::ToRay(window.tl).direction, 0); cb->frustumRayTR = Float4(Graphics3D::ToRay(window.tr).direction, 0); cb->frustumRayBL = Float4(Graphics3D::ToRay(window.bl).direction, 0); cb->size = Window::ClientRect().size; } void updateCBSky(const Float3& sunDirection) { const auto& camera = Graphics3D::GetCamera(); cbSky->cameraPos = camera.pos; cbSky->lightPos = sunDirection.normalized(); cbSky->fCameraHeight = cbSky->cameraPos.length(); cbSky->fCameraHeight2 = cbSky->fCameraHeight*cbSky->fCameraHeight; const float red = 0.680f; const float green = 0.550f; const float blue = 0.440f; const auto invWav = [](float waveLength) { return 1.0f / pow(waveLength, 4.0f); }; cbSky->invWavelength = Float3(invWav(red), invWav(green), invWav(blue)); cbSky->fInnerRadius = 1000.0f; cbSky->fInnerRadius2 = cbSky->fInnerRadius*cbSky->fInnerRadius; cbSky->fOuterRadius = 1025.0f; cbSky->fOuterRadius2 = cbSky->fOuterRadius2*cbSky->fOuterRadius2; cbSky->fScale = 1.0f / (cbSky->fOuterRadius - cbSky->fInnerRadius); const float Kr = static_cast<float>(gui.slider(L"Kr").value); const float Km = static_cast<float>(gui.slider(L"Km").value); const float Esun = 1300.0f; cbSky->fKrESun = Kr*Esun; cbSky->fKmESun = Km*Esun; cbSky->fKr4PI = Kr*4.0f*PiF; cbSky->fKm4PI = Km*4.0f*PiF; const float rayleighScaleDepth = 0.25f; cbSky->fScaleDepth = rayleighScaleDepth; cbSky->fScaleOverScaleDepth = (1.0f / (cbSky->fOuterRadius - cbSky->fInnerRadius)) / rayleighScaleDepth; cbSky->g = -0.999f; cbSky->exposure = static_cast<float>(gui.slider(L"Exposure").value); } using Uint2 = Vector2D<unsigned>; #pragma warning(push) #pragma warning(disable:4324) __declspec(align(16)) struct CB { Float4 frustumRayTL; Float4 frustumRayTR; Float4 frustumRayBL; Uint2 size; }; __declspec(align(16)) struct SkyCB { Float3 cameraPos; // The camera's current position float fCameraHeight; // The camera's current height Float3 lightPos; // The direction vector to the light source float fCameraHeight2; // fCameraHeight^2 Float3 invWavelength; // 1 / pow(wavelength, 4) for the red, green, and blue channels float fScale; // 1 / (fOuterRadius - fInnerRadius) float fOuterRadius; // The outer (atmosphere) radius float fOuterRadius2; // fOuterRadius^2 float fInnerRadius; // The inner (planetary) radius float fInnerRadius2; // fInnerRadius^2 float fKrESun; // Kr * ESun float fKmESun; // Km * ESun float fKr4PI; // Kr * 4 * PI float fKm4PI; // Km * 4 * PI float fScaleDepth; // The scale depth (i.e. the altitude at which the atmosphere's average density is found) float fScaleOverScaleDepth; // fScale / fScaleDepth float g; float exposure; }; #pragma warning(pop) PixelShader shader; ConstantBuffer<CB> cb; ConstantBuffer<SkyCB> cbSky; GUI gui; GuideArrow directionGUI; TimerMillisec directionTimer; bool guiActiveness; }; /* 参照: http://k-ichikawa.blog.enjoy.jp/etc/HP/js/sunShineAngle/ssa.html Z方向:北、X方向:東 day = [0, 365] time = [0.0, 24.0] */ Vec3 SunDirection(int day, double time, double eastLongitude, double northLatitude) { const double omega = TwoPi / 365.0; const double J = day + 0.5; const double sigma = Radians(0.33281 - 22.984*Cos(omega*J) - 0.34990*Cos(2 * omega*J) - 0.13980*Cos(3 * omega*J) + 3.7872* Sin(omega*J) + 0.03250*Sin(2 * omega*J) + 0.07187*Sin(3 * omega*J)); const double e = 0.0072*Cos(omega*J) - 0.0528*Cos(2 * omega*J) - 0.0012*Cos(3 * omega*J) - 0.1229*Sin(omega*J) - 0.1565*Sin(2 * omega*J) - 0.0041*Sin(3 * omega*J); const double T = time + (eastLongitude - 135.0) / 15.0 + e; const double t = Radians(15.0*T - 180.0); const double phi = Radians(northLatitude); const double h = Asin(Sin(phi)*Sin(sigma) + Cos(phi)*Cos(sigma)*Cos(t)); const double sinA = Cos(sigma)*Sin(t) / Cos(h); const double cosA = (Sin(h)*Sin(phi) - Sin(sigma)) / Cos(h) / Cos(phi); const double y = Sin(h); const double xz_length = Sqrt(1.0 - y*y); const double z = cosA * xz_length; const double x = sinA * xz_length; return Vec3(x, y, z); } /* 参照: http://play-siv3d.hateblo.jp/entry/jp/example/terrain */ ImageR32F CreateHeightMap() { PerlinNoise noise; ImageR32F image(512, 512); for (auto y : step(image.height)) for (auto x : step(image.width)) image[y][x].r = static_cast<float>(noise.octaveNoise(x / 256.0, y / 256.0, 8)); return image; } MeshData CreateTerrain(const ImageR32F& heightMap) { if (heightMap.isEmpty || heightMap.width != heightMap.height) { return MeshData(); } const int gridSize = heightMap.width - 1; MeshData meshData; meshData.vertices.reserve((gridSize + 1) * (gridSize + 1)); for (auto iy : step(gridSize + 1)) { for (auto ix : step(gridSize + 1)) { const double x = (ix - gridSize*0.5) / gridSize * 320; const double y = -(iy - gridSize*0.5) / gridSize * 320; const double z = heightMap[iy][ix].r * 32.0; const Vec3 position(x, 990.0 + z, y); const Vec2 uv(64.0*ix / gridSize, 64.0*iy / gridSize); MeshVertex v; v.position = position; v.texcoord = uv; meshData.vertices.push_back(v); } } auto& indices = meshData.indices; indices.resize(gridSize*gridSize * 6); for (auto y : step(gridSize)) { for (auto x : step(gridSize)) { const int index = y * (gridSize + 1) + x; indices.push_back(index + 0); indices.push_back(index + 1); indices.push_back(index + (gridSize + 1)); indices.push_back(index + (gridSize + 1)); indices.push_back(index + 1); indices.push_back(index + (gridSize + 2)); } } meshData.computeNormals(); return meshData; } /* 参照: https://gist.github.com/Reputeless/35be5b1705536cab8f43 */ class CameraManager { public: void freeCamera(double speed) { const double move = speed * (Input::KeyControl.pressed ? 20.0 : Input::KeyShift.pressed ? 5.0 : 1.0); if (Input::KeyLeft.pressed) { r -= 1.0; if (r < 0.0) r += 360.0; } if (Input::KeyRight.pressed) { r += 1.0; if (r > 360.0) r -= 360.0; } const double rad = Math::Radians(r); const double xr = Math::Sin(rad); const double zr = Math::Cos(rad); if (Input::KeyW.pressed) { state.pos.x += move * xr; state.pos.z += move * zr; } if (Input::KeyS.pressed) { state.pos.x -= move * xr; state.pos.z -= move * zr; } if (Input::KeyA.pressed) { state.pos.x -= move * zr; state.pos.z += move * xr; } if (Input::KeyD.pressed) { state.pos.x += move * zr; state.pos.z -= move * xr; } if (Input::KeyUp.pressed) { state.lookat.y += 0.5; } if (Input::KeyDown.pressed) { state.lookat.y -= 0.5; } if (Input::KeyE.pressed) { state.pos.y += move; state.lookat.y += move; } if (Input::KeyX.pressed) { state.pos.y -= move; state.lookat.y -= move; } state.lookat.x = state.pos.x + 24.0 * xr; state.lookat.z = state.pos.z + 24.0 * zr; Graphics3D::SetCamera(state); } void set(const Camera& newState) { state = newState; const Vec3 direction = (state.lookat - state.pos).normalized(); r = Degrees(Atan2(direction.z, direction.x)); Graphics3D::SetCamera(state); } private: Camera state = Graphics3D::GetCamera(); double r = 0.0; }; void Main() { Window::Resize(1280, 720); Graphics::SetBackground(Palette::Black); PS_Sky sky; CameraManager camera; camera.set(Camera({ 0, 997, 0 }, { 100, 997, 50 }, { 0, 1, 0 }, 45.0, 0.1)); const Texture textureTerrain(L"Example/grass.jpg", TextureDesc::For3D); const Mesh meshTerrain(CreateTerrain(CreateHeightMap())); Graphics3D::SetAmbientLight(ColorF(0.2, 0.3, 0.4)); Graphics3D::SetFog(Fog::SquaredExponential(ColorF(0.5, 0.6, 0.8), 0.002)); //GuideArrow用 Graphics3D::SetLightForward(Light().Directional(Float3(0, 1, 0).normalized())); Graphics3D::SetAmbientLightForward(ColorF(0.3)); Window::SetTitle(L"大気散乱のシミュレーション | [H]キーでGUIの表示/非表示切り替え"); while (System::Update()) { camera.freeCamera(0.5); if (sky.pushedResetButton()) { camera.set(Camera({ 1, 997, 1 }, { 1 + 100, 997, 1 }, { 0, 1, 0 }, 45.0, 0.1)); sky.reset(); } sky.update(); Graphics3D::SetLight(0, Light().Directional(sky.sunDirection())); meshTerrain.draw(textureTerrain); } }
参考文献
(1) GPU Gems 2 日本語版、p233 - p246
(2) http://hexadrive.sblo.jp/article/47534333.html、ヘキサドライブ日記 大気散乱シミュレーション -Atmospheric Scattering-
(3) http://k-ichikawa.blog.enjoy.jp/etc/HP/js/sunShineAngle/ssa.html、音と色と数の散歩道、太陽の高度と方位角(JavaScriptによる簡易版)
(4) http://play-siv3d.hateblo.jp/entry/jp/example/terrain、Play Siv3D!、サンプル | 地形
(5) https://gist.github.com/Reputeless/35be5b1705536cab8f43、GitHub Gist、FreeCanera.cpp