레이블이 Graphic인 게시물을 표시합니다. 모든 게시물 표시
레이블이 Graphic인 게시물을 표시합니다. 모든 게시물 표시

2021년 4월 7일 수요일

Sky Atmospheric Scattering

참고한자료

 https://developer.nvidia.com/gpugems/gpugems2/part-ii-shading-lighting-and-shadows/chapter-16-accurate-atmospheric-scattering

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 대용으로 쓰면 됨.




 아 잘된다.




List