참고한자료
https://hal.inria.fr/inria-00288758/document
https://github.com/OzgurCerlet/precomputed_atmospheric_scattering_dx11
https://old.cescg.org/CESCG-2009/papers/PragueCUNI-Elek-Oskar09.pdf
https://github.com/Scrawk/Brunetons-Improved-Atmospheric-Scattering
내가 따로 주석단 거
https://github.com/Mona04/my-try-to-understand-Brunetons
1. Sean O'Neil
단일 Scattering 만을 고려한 런타임 중에 가능한 기초적인 구현이 되어있음.
Precompute 된 데이터를 Scale 함수로 근사해 들고옴.
//GPUGem2 //https://developer.nvidia.com/gpugems/gpugems2/part-ii-shading-lighting-and-shadows/chapter-16-accurate-atmospheric-scattering // Calculates the Mie phase function float getMiePhase(float fCos, float fCos2, float g, float g2) // g = (-1, 1) { return 1.5 * ((1.0 - g2) / (2.0 + g2)) * (1.0 + fCos2) / pow(1.0 + g2 - 2.0 * g * fCos, 1.5); } // Calculates the Rayleigh phase function float getRayleighPhase(float fCos2) // phase with g = 0 { return 0.75 + 0.75 * fCos2; } // The number of sample points taken along the ray static const int nSamples = 2; static const float fSamples = (float) nSamples; // The scale depth (the altitude at which the average atmospheric density is found) static const float Kr4PI = 0.0025f * 4.0f * pi; static const float Km4PI = 0.0010f * 4.0f * pi; static const float KrESun = 0.0025f * 10; // sun = 20 static const float KmESun = 0.0010f * 10; static const float g = -0.99f; static const float g2 = g * g; static const float ScaleDepth = 0.25; static const float InvScaleDepth = 1.0 / ScaleDepth; static const float InnerRadius = 10.0f; static const float OuterRadius = 10.25f; static const float HeightScale = 1 / (OuterRadius - InnerRadius); static const float ScaleOverScaleDepth = HeightScale / ScaleDepth; static const float3 waveLength = { 0.650f, 0.570f, 0.475f }; // 0.001 = 1nm static const float3 waveLength4_inv = { 1/pow(waveLength[0], 4), 1/pow(waveLength[1], 4), 1/pow(waveLength[2], 4) }; // 0.001 = 1nm // The scale equation calculated by Vernier's Graphical Analysis // ScaleDepth, OuterRadius - InnerRadius 가 0.25 로 일치할 때만 기능함. float scale(float fCos) { float x = 1.0 - (fCos < 0 ? -fCos : fCos); return ScaleDepth * exp(-0.00287 + x * (0.459 + x * (3.83 + x * (-6.80 + x * 5.25)))); } void Scatter_Ground(float3 v_pos, float3 sun_dir, out float3 rayleigh, out float3 mie) { // Get the ray from the camera to the vertex, and its length (which is the far point of the ray passing through the atmosphere) float3 v_ray = v_pos; float fFar = length(v_ray) * 0.25; // Calculate the ray's starting position, then calculate its scattering offset float3 ray_start = float3(0, InnerRadius + 0.001, 0); float fHeight = length(ray_start); float fDepth = exp(ScaleOverScaleDepth * (InnerRadius - fHeight)); // -h / h0 float fStartAngle = dot(v_ray, ray_start) / fHeight; // 땅을 보면 0이고 하늘을 보면 1이 되게 하려고 float fStartOffset = fDepth * scale(fStartAngle); // Initialize the scattering loop variables float fSampleLength = fFar / fSamples; float fScaledLength = fSampleLength * HeightScale; float3 sampleRay = v_ray * fSampleLength; float3 samplePoint = ray_start + sampleRay * 0.5f; // Now loop through the sample rays float3 front_color = float3(0.0, 0.0, 0.0); for (int i = 0; i < nSamples; i++) { float fHeight = length(samplePoint); float fDepth = exp(ScaleOverScaleDepth * (InnerRadius - fHeight)); float fLightAngle = dot(-sun_dir, samplePoint) / fHeight; float fCameraAngle = dot(v_ray, samplePoint) / fHeight; float fScatter = (fStartOffset + fDepth * (scale(fLightAngle) - scale(fCameraAngle))); float3 attenuate = exp(-fScatter * (waveLength4_inv * Kr4PI + Km4PI)); // rayleigh 만 wave 에 영향을 받음 front_color += attenuate * (fDepth * fScaledLength); samplePoint += sampleRay; } // Finally, scale the Mie and Rayleigh colors and set up the varying variables for the pixel shader rayleigh = front_color * (waveLength4_inv * KrESun); mie = front_color * KmESun; }
float4 PS_Dome(PixelPositionColorUV input) : SV_Target { float3 sun_dir = normalize(sun_pos); float3 v_dir = normalize(input.pos.xyz); float3 rayleigh, mie; Scatter_Ground(v_dir, -sun_dir, rayleigh, mie); float fCos = dot(-sun_dir, v_dir); float fCos2 = fCos * fCos; float3 color = rayleigh + getMiePhase(fCos, fCos2, g, g2) * mie; return float4(color, 1); }
2. Brunetons
위 과정을 따라가게 됨.
근데 복잡하니까 Single Scattering 까지는 다른 논문의 수식으로 설명하겠음.
선이 두줄로 된 글자가 우리가 텍스쳐로 저장할 데이터고 검은색으로 꽉찬 대문자가 함수임.
T 는 Transmittance, E 는 Irradiance, S 는 Scattering 으로 우리는 S 를 구하는걸 목표로함.
1. Transmittance(Optical Length)
1.Parameter
우리는 2D Texture 을 만들어야하고 우리가 필요한 인자는 view-zenith 의 각도와 현재의 높이임. 우리가 지구를 완전한 구로 가정하므로 z 축의 값은 필요없음.
논문을 보면 view-zenith 각도를 인자로 직접 쓸 때 우리가 바라보는 대상의 거리에 따른 인자의 변화가 극명하게 나타나지 않는다고 함.(델타가 지평선 부분에서 너무 작음) 그래서 이제 설명한 것처럼 인자를 바꿔서 ( 거리를 인자로 바로 씀 ) 이를 해결함.
우선 각 변수를 소개함.
top 은 지구 중심과 대기끝까지의 거리, bottom 은 지구 중심과 대기 시작과의 거리로 고정된 값임. r 는 현재 관측자의 높이, mu 는 zenith 와 빛이 들어오는 방향(보통 vertex)에 대한 cos 으로 이 둘은 Input 이 될 중요한 인자임. rho 는 관측자의 지평선이 땅과 닿는 지점과의 거리. H 는 sqrt(top^2 - botton^2) 으로 mu, r 을 구하는데 쓰일 값임.
d_min 은 대기 끝과 관측자의 최단거리로 top - r 임. d_max 는 최대거리로 지평선이 대기 끝과 닿는 곳까지의 거리인 H + dh0 임. d 는 view vector 가 대기끝에서부터 관측자까지로 잘린 거리로 위 그림의 주황색 호 사이에서 가능함.
uv 가 인자로 들어오는데 uv_x 값은 d_min 과 d_max 사이를 보간하는 인자가되고 uv_y 는 H 값을 0~1 사이로 쪼갠 값이 되어 H * uv_y = rho 가 됨. d_min < d < d_max 는 위 그림상 자명하므로 uv_x 가 d 를 표현할 수 있는건 쉽게 알 수 있음. uv_y 에선 rh0 가 H 와 같으면 r 은 대기끝으로 가야만하고 0 이 되면 지표면에 있게 되므로 충분히 관측자의 높이를 알 수 있게 됨.
cosθ 는 코사인 제2법칙을 사용해서 구하는데 앞에서 설명했듯이 180-θ 부분을 구하고 cos(180-θ) = -cosθ 이므로 여기에 음수처리를 하게 됨. 코드에선 다음과 같음.
(H^2 - rho^2 - d^2)/(2*r*d)
= -(rho^2 + d^2 - H^2) / (2 * r * d)
= -(rho^2 + d^2 - top^2 + bottom^2) / (2 * r * d)
= -(r^2 + d^2 - top^2) / (2 * r * d)
즉 d, r, top 이 이루는 삼각형에 대해서 d, r 사이의 코사인 값을 얻고 이를 음수처리한 것.
mu, r 로 d 를 구해야할 경우가 자주 있어 이것도 설명함. 이번의 d 는 두가지 종류가 있음. 하나는 대기 끝에 닿는 경우이고 하나는 대지에 닿는 경우임. 전자는 하늘을 보았을 때고 후자는 땅에 반사된 빛을 계산할때 등에 쓰임. 위에선 d 가 전자의 의미였음에 참고.
이를 구하기 위해선 원점은 지구의 중심이고, x 축이 view vector 이고(즉 d 를 이루는 벡터), y 축은 그것과 수직이며 양의 방향이 관측자 쪽이라고 생각해야함. 이 좌표축에서 예를들어 (d, 0) 은 원점에서 위 그림의 d 가 가진 벡터의 방향만큼 이동항 것이 됨. 그럼 (r * mu, r * sqrt(1-mu^2)) 의 위치는 관측자의 위치가 될 것임. 여기서 (d + r * mu, r * sqrt(1-mu^2) ) 는 대지 똑은 대기끝과의 접점이 될 것임.
대기와의 접점일 때는 d 를 위 식에서부터 유도할 수 있음.
top^2 = r^2(1-mu^2) + (d + r*mu)^2
sqrt(top^2 + r^2(mu^2-1)) = sqrt(discriminant) = d + r*mu
대지와의 접점일 때는 d + r*mu 가 음수이므로 sqrt 할 때 - 를 붙여야해서 다음과 같음.
sqrt(top^2 + r^2(mu^2-1)) = sqrt(discriminant) = -(d + r*mu)
2. Compute
Lambert-Beer law 에 따른 공식임. Rayleigh, Mie, 기타 흡수 등을 각각 구해서 더해줘야한다. 베타 함수는 Scatter 상수로 Rayleigh 의 경우 파란색같은 고파장이 더 크게 되어 하늘이 파란 이유가 되기도 한다. 로 함수는 산란 종류에 따른 입자의 농도로 지표면으로부터의 높이에 따른 식이다.
이게 베타에 대한 식인데 N, n 모두 해수면 기준으로 한 값이다. 람다는 광선의 파장 길이로 Rayleigh 위에서 말한 특징이다. 이 식들은 모두 상수이므로 Brunetones 은 구현할 때 상수로 빼 두었고 논문에서 베타는 알파 베타 로 모두 합쳐져 하나로 표현된다.
2. Direct Irradiance
Irradiance 는 두가지로 나뉨. 에너지원으로부터 바로오는 빛에 대한 Direct Irradiance 와 어떤 물체와 부딪힌 후에 오는 경우인 Indirect Irradiance 가 그것. 이 둘이 구하는 것이 다름. 전자는 아니지만 후자의 경우 지표면으로부터의 높이와 노멀값이 모두 필요하게 됨. 하지만 계산 편의상 부딪히는 면이 평평하다고 치고 높이만 가지고 값을 구하게 됨.
Indirect Irradiance 는 계산해야할게 많지만 Direct Irradiance 는 위에서 구한 Transmittance 만으로 구할 수 있으므로 이걸 먼저 계산하게 됨. Multi Scattering 에서 쓰는 값임.
정확히는 위 식을 따르는데, 땅인 x0 에서 반사되어서 x 지점의 빛의 양을 R 이라고 하게 됨. x 는 구할때마다 달라지고 x0 은 지금 바로 알 수 있음. 그래서 아래 크시 함수부분만 여기서 구함. 그리고 지금은 태양광에서 빛이 바로오는 거만 고려하므로 적분은 안해도 됨.
float3 ComputeDirectIrradiance(Texture2D transmittance_texture, float r, float mu_s) { // 우리 시야를 구로 칠 때 태양이 차지하는 범위의 각도의 sin 혹은 지평선에서 태양 지름위까지의 cos
float alpha_s = _sun_angular_radius; // cos_factor 는 0 ~ 1 사이의 값이 됨. // 지평선 너머로 갈 시 태양이 지표에 가리므로 이걸 구현하기 위해서 // 0 < mu_s < alpha_s 의 경우 0~alpha_s 의 값이되 0으로 갈 때 급격히 없어지게 구현함. float average_cosine_factor = (mu_s < -alpha_s) ? 0.0 : ((mu_s > alpha_s) ? mu_s : (mu_s + alpha_s) * (mu_s + alpha_s) / (4.0 * alpha_s)); // sun solid angle 이 작으므로 transmittance 를 상수로 취급함. ( 원래는 적분해야함 ) return _solar_irradiance * average_cosine_factor * GetTransmittanceToTopAtmosphereBoundary(transmittance_texture, r, mu_s); }
sun_angular_radius 는 mu 를 구한 방식으로 치면 태양 반의 시야각만큼 지평선에서 올라온 cos 을 말함. 혹은 태양 반의 시야각만큼의 sin 을 말함. cos(90-x) = sin(x) 이므로 이 둘은 당연히 같음.
우리가 친숙한 Lambertian surface 모델에서 Diffuse 구하는 거와 같음. 광원이 나에게 오는 각도에 따라 수직일 땐 온전히 수평에 가까울 땐 0에 가깝되 되는 - 정확히는 cos(각도) - 그것임. 우리는 평평하다고 가정을 했으므로 dot( (0, 1, 0), 태양 각도) = cos(태양각도) = mu 가 됨 (mu 는 parametar 참고).
3. Single Scattering
대기 중에 단 한번의 산란을 한걸 계산하는 것. 그래서 카메라 위치(p)와 대기 중의 산란 지점(q) 그리고 q 와 태양 간의 벡터와 대기 끝과의 접점(i?), 이 세 점 사이를 지나는 빛에 대해서 scattering 을 계산하게 됨.
1. Parameter
float mu_s = dot(camera, sun_direction) / r; float nu = dot(view_ray, sun_direction);
총 4개의 인풋이 필요하며 3D Texture 가 공간의 한계라 한 축으로 두개의 인자에 할당해서 사용함. 그 인자는 mu, r 에 더해서 mu_s, nu 가 추가됨. mu_s 와 nu 는 위 코드대로의 의미인데 둘이 합쳐서 view-run 사이의 각도를 확정함. mu, r 은 transmittance 와 다르게 ray 가 하늘과 닿을 때랑 지면과 닿을 때를 모두 포괄하게 됨.
void GetRMuMuSNuFromScatteringTextureUvwz(const in float4 uvwz, out float r, out float mu, out float mu_s, out float nu, out bool ray_r_mu_intersects_ground) { float H = sqrt(_top_radius2 - _bottom_radius2); float rho = H * GetUnitRangeFromTextureCoord(uvwz.w, SCATTERING_TEXTURE_R_SIZE); r = sqrt(rho * rho + _bottom_radius2); if (uvwz.z < 0.5) { float d_min = r - _bottom_radius; float d_max = rho; // 0<= z <0.5 이므로 d_max ~ d_min 의 값이 됨. // 각도가 천정으로부터 땅으로 내려오는걸 생각하면 땅에 부딪히는 지점부터 발끝을 볼때까지가 되는 것 float d = d_min + (d_max - d_min) * GetUnitRangeFromTextureCoord(1.0 - 2.0 * uvwz.z, SCATTERING_TEXTURE_MU_SIZE / 2); mu = d == 0.0 ? float(-1.0) : ClampCosine(-(rho * rho + d * d) / (2.0 * r * d)); // (bottom2 - r2 - rho2) / 2rd = (-d2 -rh02) / 2rd ray_r_mu_intersects_ground = true; } else { float d_min = _top_radius - r; float d_max = rho + H; float d = d_min + (d_max - d_min) * GetUnitRangeFromTextureCoord(2.0 * uvwz.z - 1.0, SCATTERING_TEXTURE_MU_SIZE / 2); mu = d == 0.0 ? float(1.0) : ClampCosine((H * H - rho * rho - d * d) / (2.0 * r * d)); ray_r_mu_intersects_ground = false; } float x_mu_s = GetUnitRangeFromTextureCoord(uvwz.y, SCATTERING_TEXTURE_MU_S_SIZE); float d_min = _top_radius - _bottom_radius; float d_max = H; float A = -2.0 * _mu_s_min * _bottom_radius / (d_max - d_min); float a = (A - x_mu_s * A) / (1.0 + x_mu_s * A); float d = d_min + min(a, A) * (d_max - d_min); // A = _bottom / (d_max - d_min) // 0 => d_min + _bottom_radius = top_radius // 1 => d_min mu_s = d == 0.0 ? float(1.0) : ClampCosine((H * H - d * d) / (2.0 * _bottom_radius * d)); // top2 - bottom2 - d2 / 2dbottom = -(d2 + bottom2 - top2)/2dbottom nu = ClampCosine(uvwz.x * 2.0 - 1.0); }
r, mu 는 분기로 나뉜걸 빼면 특별한게 없음. 코사인 제 2법칙으로 각도를 알아내게 됨.
mu_s, nu 가 문제인데 mu_s 의 의미는 위의 설명 대로지만 uv 에서 얻는 과정은 정확히 어떤 과정인지 모르겠음. 위 코드대로 계산하면 다음과 같이 나옴.
x_mu_s = 1 => d = d_min + bottom
x_mu_s = 0.87 = 1 - (875.67 - 60)/6360 = (1-(d_max-d_min)(bottom)) / 2 => d = d_max
x_mu_s = 0 => d_min
그래서 x_mu_s 는 0 일 때 cos(0), 0.87 일 대 cos(90), 1 일때 cos(90+xxx) 가 되고 따라서 인자의 정의역인 (0~1) 의 대부분을 각도과 0 ~ 90 사이에 배치해 둠을 알 수 있음. mu_s_min 이 기본으로 -0.5 로 되어있는데 이게 만약 -1 이면 cos(90+xxx) 가 cos(180) 이 될거임.
정리하면 mu_s_min 에 따라서 d 가 d_max 인 H 를 넘어서 빨간 선을 따라 증가하게 되고, 이 초과범위는 0.1x 정도의 도메인의 영역을 차지함.
그런데 그냥 대충 식을 끼워맞춰서 저런 식이 나왔는지 어쩐지는 모르겠음.
nu 는 mu_s 가 cos 값 즉 우함수라 빛의 방향을 온전히 지정할 수 없어서 사용됨.
nu = clamp(nu, mu * mu_s - sqrt((1.0 - mu * mu) * (1.0 - mu_s * mu_s)), mu * mu_s + sqrt((1.0 - mu * mu) * (1.0 - mu_s * mu_s)));
이게 위함수 뒤에 적용이 되는데 자세한건 및에서 설명함.
2. 산란 지점 구하기
우리는 위에서 대기중의 점 p 에서 대기 끝까지의 Transmittance 를 구할 수 있었음. 이걸 이용해서 우리는 대기중의 점 p 와 대기 중의 점 q 사이의 Transmittance 도 구할 수 있음.
이를 구하기 위해서 r, mu, d, mu_s, nu 가 인풋으로 들어옴. transmittance 에서는 d 가 당연히 p 와 대기 끝의 거리였지만 이번엔 임의의 길이가 됨. 이 길이로 우리는 q 의 위치를 특정할 수 있음. 그리고 바로 그 특정하기 위한 값인 mu_s_d 와 r_d 를 구해야함.
r_d 는 코사인 2법칙을 쓰면 간단하게 구해짐. r 과 d 그리고 mu 를 알고있기 때문임. 이때 r 과 d 사이의 코사인 값은 -mu 임에 주의.
mu_s_d 는 복잡함. 우선 위 그림을 보면 우리가 원하는건 하늘색 직선(q)과 노란색 직선(태양, s) 이 이루는 cos 값임. 그리고 우리가 갖고 있는 mu_s 는 태양과 p 가 이루는 cos 값임. 그래서 우리는 (r*mu_s + nu_s*d)/r_d 를 해서 이를 해결함. 이는 가장 아래쪽의 노란 선을 기준으로 코사인을 하면 이해가 쉬움. 이때 위 그림에선 nu_s * d 가 음수임에 주의.
여기서 주의해야할 것은 nu_s 가 view_ray 와 태양빛의 벡터 사이의 각도에 따라서 나타내는 각도가 다르다는 점임. 위 그림에선 mu = cos(theta_d) 라고 하면, nu_s = cos(theta_d + theta_s) 이어야 함. 하지만 빛의 방향이 반대가 된다면 nu_s = cos(theta_d - theta_s) 여야만 함. 직접 위 그림에서 빛만 / 방향으로 그려보면 왜 그런지 답이 나옴.
그래서 nu 에 clamp 를 거는 이유 중 하나를 할 수 있음. 그 식을 풀면
nu = clamp(nu, cos(theta_d + theta_s), cos(theta_d - theta_s)) 이기 때문.
이 clamp 사이의 값은 지금까지 2D 차원에서 보고 있었기 때문에 설명하지 않았음. 만약 빛이 3D 차원이라면 보라색 선인 view_dir 과 다른 각도에서 빛을 그릴 수 있음. 이는 cos(camera, sun_dir) 인 mu_s 에선 상관이 없고 cos(view_dir, sun_dir) 인 nu 에서 관계가 있음. mu_s 가 고정되었다고 했을 때 nu 는 최대 위 그림처럼되고 최소론 빛이 / 방향으로 대체되었을 때가 될 거임. 그리고 빛이 위 그림의 2D 차원이 아닌 다른 축으로 기울어져 있다면 nu 의 값은 theta_d +theta_s ~ theta_d - theta_s 사이가 될 것임.
3. 두 점 사이의 Transmittance
float3 GetTransmittance(const in Texture2D transmittance_texture, float r, float mu, float d, bool ray_r_mu_intersects_ground) { // mu 는 d 와 r 이 r_d 가 이루는 각에 대한 -cos 인걸 참고 float r_d = ClampRadius(sqrt(d * d + 2.0 * r * mu * d + r * r)); float mu_d = ClampCosine((r * mu_s + d * nu) / r_d); if (ray_r_mu_intersects_ground) { return min(GetTransmittanceToTopAtmosphereBoundary(transmittance_texture, r_d, -mu_d) / GetTransmittanceToTopAtmosphereBoundary(transmittance_texture, r, -mu), (float3) (1.0)); } else { return min(GetTransmittanceToTopAtmosphereBoundary(transmittance_texture, r, mu) / GetTransmittanceToTopAtmosphereBoundary(transmittance_texture, r_d, mu_d), (float3) (1.0)); } }
mu_s_d, r_d 를 알았으니 각 점에서 대기 끝까지의 Transmittance 의 차를 구하면 그 두 점 사이의 값이 나올 것임. 이때 우리는 위 코드처럼 값을 나누게 되는데 이건 이 값이 지수에 들어가기 때문임.
위 식은 Single Scattering 을 구하는 식임. 지금 주의해서 볼 곳은 exp 안쪽임. 위 식에서 Pa 는 관측자, P 는 산란 지점, Pc 는 태양광이 대기끝과 접하는 지점으로 Pa~P, P~Pc 의 Transmittance 를 더하고 있음 ( exp(-transmittance) 임을 참고). 그런데 우리가 텍스쳐에 저장한 값인 exp(-density) 이므로 더하는게 아니라 나누는 걸로 구해야함. 우리는 위 식 표현을 빌리면 Pa~Pc 와 P~Pc 사이의 Transmittance 를 구해서 나누므로써 같은 식을 구하게 됨.
덧붙여 위 코드에서 인자로 대지와 충돌한지 받는 이유는 r, mu 의 resolution 이 낮아서 지평선 근처에선 오작동을 하고, caller 가 이미 그 값을 간단히 알 고 있기 때문임.
4. 임의의 점과 태양 사이의 Transmittance
원래는 태양의 지름에 따라서 측점 지점과 태양 사이는 원뿔의 공간이 있기 때문에 그 내부의 transmittance 를 적분해야 할 것임. 하지만 태양이 멀리있으므로 광선은 평행하다고 해도 무관함. 그래서 Transmittance 는 태양의 크기에 대해서 상수로 취급함. 그리고 태양이 지표면에 가까이가서 가려질 때만 이를 처리하게 됨.
float3 GetTransmittanceToSun(const in Texture2D transmittance_texture, float r, float mu_s) { //cos(theta_h + alpha) <= mu <= cos(theta_h - alpha) //cos(theta_h)cos(alpha) - sin(theta_h)sin(alpha) <= mu <= cos(theta_h)cos(alpha) + sin(theta_h)sin(alpha) //-sin(theta_h)sin(alpha)<= mu - cos(theta_s)cos(alpha) <= +sin(theta_h)sin(alpha) // cos(alpha) 는 1에 가까우니까 1로 침 float sin_theta_h = _bottom_radius / r; // sin(180-theta_h) = sin(theta_h) float cos_theta_h = -sqrt(max(1.0 - sin_theta_h * sin_theta_h, 0.0)); return GetTransmittanceToTopAtmosphereBoundary(transmittance_texture, r, mu_s) * smoothstep(-sin_theta_h * _sun_angular_radius, sin_theta_h * _sun_angular_radius, mu_s - cos_theta_h); }
theta_h 는 지평선까지의 각도이고 태양과의 각도는 theta_s 로 표현함. 즉 mu_s = cos(theta_s) 임. 그리고 sun_angular_raidus = sin(alpha) 는 시야에서 태양이 차지하는 각도에 대한 사인 값임. 그러면 다음과 같은 부등식 사이에서 태양이 부분적으로 존재함.
cos(theta_h + alpha) <= mu_s <= cos(theta_h - alpha)
cos(theta_h)cos(alpha) - sin(theta_h)sin(alpha) <= mu_s <= cos(theta_h)cos(alpha) + sin(theta_h)sin(alpha) // cos(a+-b) 를 풀어쓴것
-sin(theta_h)sin(alpha)<= mu_s - cos(theta_s)cos(alpha) <= +sin(theta_h)sin(alpha)
-sin(theta_h)sin(alpha)<= mu_s - cos(theta_s) <= +sin(theta_h)sin(alpha)
// cos(alpha) 는 1에 가까우니까 1로 침
5. Precompute
이제 위 함수를 적용하면 됨.
베타 함수는 Scatter 상수이고 로 함수는 고도에 따른 기체의 농도로 위에서 설명함.
I 는 빛의 강도임.
위 수식에선 Pa~Pc 사이의 Transmittance 를 P~Pc, Pa~P 사이의 Transmittance 를 더하는데, Pa~P 가 두 점 사이의 Transmittance 가 될 것이고 P~Pc 가 임의의 점과 태양 사이의 Transmittance 가 될 것임.
F 는 Phase 함수로 세타 방향으로 산란되었다고 했을 때 그 정도를 나타내는 근사 함수임. g 는 비대칭 정도를 말하는데 Rayleigh 의 경우 g = 0 이고 Mie 의 경우 Brunetons 는 0.8 을 사용함. 그리고 코드 내에서 보면 여기에 1/4pi 를 같이 넣어놨으니 주의.
그리고 Phase 는 GetScatter 에서 계산되고 Mie, Rayleigh 텍스쳐에 저장할 땐 빼고 넣음. 왜냐면 각 각도에 맞춘 데이터는 해상도 등으로 실제 각도에 부드럽게 적용되기 힘들기 때문.
4. Multi Scattering
이번엔 한번의 산란이 아니라 여러번의 산란을 거쳐서 한 점에 도달한 빛의 양을 구함. 이때 마지막 산란이 땅에 의한 산란인 경우는 배제함. 왜냐하면 이건 렌더링 때 알베도 등을 고려해야하기 때문임. 그 이전에 땅과 산란한 경우 편의 상 값을 넣은 Uniform Albedo 를 써서 값을 구하게 됨.
그리고 계산에 사용되는 빛인 L 은 Single Scattering 에서는 태양광을 그대로 썼지만 2번 이상의 산란에선 그렇지 않음. 그리고 이는 위의 식을 따르게 됨. 이때 R 는 지표 등에 반사된 빛이고 S 는 Scattering 을 말함.
위식은 light_shaft 때문에 있던거 같은데 그건 고려 안했음.
1. Density
앞에서 Single Scattering 을 구했는데 이번엔 한 점을 향에 모든 방향에서 오는 빛을 합쳐야하기 때문에 위 식은 너무 복잡함.
그래서 exp(-t) 부분과 나머지 부분으로 적분을 나눠버림.
바로 이처럼 말임.
여기서 J 가 지금 구할 Density 가 됨.
이때 L 은 땅에서의 반사인 R + 지금까지의 Scattering 합 이 될거임. 이때 R 에 관해서, 지금까지 구한건 크시 함수 부분이므로 나머지 부분(Transmittance, 1/pi) 을 여기서 곱해주게 됨.
float3 omega = float3(sqrt(1.0 - mu * mu), 0.0, mu); //view_direction //dot(view, sun) - dot(zenith, sun) * dot(zenith, view) / omega.x //= (dot(view, sun) - sun_z * view_z) / view_x ... view 는 가정상 y 값이 없음. //= (view_x * sun_x + view_z * sun_z - sun_z * view_z) / (view_x) //= sun_x
참고로 view_direction 에 해당하는 omega 는 위처럼 구해지는 것.
2. Indirect Irradiance
Direct Irradiance 구하는 것과 똑같이 크시함수만 여기서 구함. 이젠 모든 방향에서의 빛을 구해야해서 2pi 구간에서 적분하게 됨.
이때 L 은 기존의 R 은 영향이 적다고 판단되어 Scattering 만 가져와서 구하게 됨.
3. Multi Scattering
J , T는 모두 계산을 해놨기 때문에 적분하면서 필요한 값 읽으면 됨.
5. Rendering
위에서 구한 scattering, single_mie_scattering, irradiance, transmittance 로 나머지를 구할 수 있음. 구현한거 보면 scattering 에는 rayleigh 만 여러번 적용이 되어 있어서 scattering 은 multi_rayleigh_scattering 같은 느낌임.
예제 코드의 ps 는 구와 하늘와 태양과 땅까지 한번에 다 합쳐서 lerp 로 해결함.
그 부분을 따로 떼서 적용하면 쉽게 적용 가능함.
1. Sky
float r = length(camera); float rmu = dot(camera, view_ray); // rmu = r*mu // rmu*rmu - r*r + top2 = discriminant // r^2(1-mu^2) + (d+rmu)^2 = top^2 // camera 가 우주에 있으면 d+rmu < 0 일 것이므로 d+rmu = -discriminant float distance_to_top_atmosphere_boundary = -rmu - sqrt(rmu * rmu - r * r + _top_radius2); if (distance_to_top_atmosphere_boundary > 0.0) { camera = camera + view_ray * distance_to_top_atmosphere_boundary; r = _top_radius; rmu += distance_to_top_atmosphere_boundary; } else if (r > _top_radius) { transmittance = (float3) (1.0); return (float3) (0.0); }
우주 공간을 고려해서 위와 같은 보정이 들어가게됨.
2. Ground
InScatter 랑 Irradiance 를 갖게 되는데 전자는 멀리 있는 물체가 뿌옇게 나오는 이유가 됨.
Irradiance 는 Diffuse 대용으로 쓰면 됨.
아 잘된다.