rayMarching renderer population

rayMarching_population.glsl

//This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
//Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
//associated documentation files (the "Software"), to deal in the Software without restriction,
//including without limitation the rights to use, copy, modify, merge, publish, distribute,
//sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is
//furnished to do so, subject to the following conditions:

//The above copyright notice and this permission notice shall be included in all copies or
//substantial portions of the Software.

//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
//INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
//PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
//FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
//OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
//IN THE SOFTWARE.

//---------------------------------------------
//---------- renderer.vert
//---------------------------------------------
#version 460 core

uniform vec4 vertices[] = vec4[](vec4(-1.0, -1.0, 0.5, 1.0),
                                   vec4( 1.0, -1.0, 0.5, 1.0),
                                   vec4(-1.0,  1.0, 0.5, 1.0),
                                   vec4( 1.0,  1.0, 0.5, 1.0));

void main(void)
{
    gl_Position = vertices[gl_VertexID];
}


//---------------------------------------------
//---------- renderer.frag
//---------------------------------------------

//based on https://www.shadertoy.com/view/3lsSzf

#version 460 core

layout (location = 0) out vec4 fragColor;
layout (depth_any) out float gl_FragDepth;

uniform vec2 iResolution;
uniform mat4 pmvInversed;
uniform mat4 pmv;
uniform mat4 modelview;
uniform mat4 projection;
uniform vec3 sunlight;
uniform int gridsize;
uniform float scaleZ;
uniform int indexMax;
uniform vec3 materialColor1;// = vec3(0.74, 0.81, 0.330);
uniform vec3 materialColor2;// = vec3(0.927, 0.3, 0.1);
uniform vec3 materialColor3;// = vec3(0.74, 0.81, 0.330);
uniform vec3 materialColor4;// = vec3(0.927, 0.3, 0.1);
uniform float time;
uniform float marchingUnit;

struct GridPopu { //8 Bytes 위치는 gl_vertexID 로 찾는다.
	uint sgg; //4	
	float popu; //4
};

layout (binding = 0) uniform sampler2D sColor;
layout (binding = 1) uniform sampler2D sDepth;

layout (binding = 2) buffer popuBuffer {
	GridPopu popu[];
};

vec4 groundRendered;
float worldScale =1.0f; 
//지구가 보이는 우주의 경우 100.0;
//-10~ 10.0의 작은 공간에서 0.001;

float iTime = 0.0f;
//------------------------------------------------------------------

float hash31(vec3 p){
    return fract(sin(dot(p, vec3(12.989, 78.233, 57.263)))*43758.5453);
}

float sdBox( vec3 p, vec3 b )
{
  vec3 q = abs(p) - b;
  return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}

// http://iquilezles.org/www/articles/smin/smin.htm
float smin( float a, float b, float k )
{
    float h = max(k-abs(a-b),0.0);
    return min(a, b) - h*h*0.25/k;
}

// http://iquilezles.org/www/articles/smin/smin.htm
vec2 smin( vec2 a, vec2 b, float k )
{
    float h = clamp( 0.5+0.5*(b.x-a.x)/k, 0.0, 1.0 );
    return mix( b, a, h ) - k*h*(1.0-h);
}

// http://iquilezles.org/www/articles/smin/smin.htm
float smax( float a, float b, float k )
{
    float h = max(k-abs(a-b),0.0);
    return max(a, b) + h*h*0.25/k;
}

// http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
float sdSphere( vec3 p, float s )
{
    return length(p)-s;
}

// http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
float sdEllipsoid( in vec3 p, in vec3 r ) // approximated
{
    float k0 = length(p/r);
    float k1 = length(p/(r*r));
    return k0*(k0-1.0)/k1;
}

float sdVerticalCapsule( vec3 p, float h, float r )
{
  p.z -= clamp( p.z, 0.0, h );
  return length( p ) - r;
}

vec2 sdStick(vec3 p, vec3 a, vec3 b, float r1, float r2) // approximated
{
    vec3 pa = p-a, ba = b-a;
	float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
	return vec2( length( pa - ba*h ) - mix(r1,r2,h*h*(3.0-2.0*h)), h );
}

// http://iquilezles.org/www/articles/smin/smin.htm
vec4 opU( vec4 d1, vec4 d2 )
{
	return (d1.x<d2.x) ? d1 : d2;
}


float opCheapBend(  vec3 p, float innerR, float outerR )
{
    const float k = 0.00000001; // or some other amount
    float c = cos(k*p.x);
    float s = sin(k*p.x);
    mat2  m = mat2(c,-s,s,c);
    vec3  q = vec3(p.x, m*p.yz);

    float z = dot(q, vec3(0,0,1));
    
    float xy = abs(length(q.xy - vec2(0,0)) - (innerR + outerR)/2.0f) - (outerR-innerR)/2;
    if (xy<0) xy = 0.0f;
    return sqrt((z*z)+(xy*xy));   
}


float sdCappedCylinder( vec3 p, float h, float r )
{
  vec2 d = abs(vec2(length(p.xy),p.z)) - vec2(r,h);
  return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}


float opUnion( float d1, float d2 ) { return  min(d1,d2); }

float opSubtraction( float d1, float d2 ) { return max(-d1,d2); }

float opIntersection( float d1, float d2 ) { return max(d1,d2); }


#define BASEX 1074000
#define MAXX 1116000
#define BASEY 1728000
#define MAXY 1785000

uvec2 basexy = uvec2(1074000,1728000); //~(1116000, 1785000)
uint gridMaxX = (1116000-basexy.x+(gridsize-1))/gridsize;


vec2 map(vec3 pos)
{
    int x0 = (int(pos.x) - int(basexy.x))/gridsize;
    int y0 = (int(pos.y) - int(basexy.y))/gridsize;
   
    float bar0 = 1000000000.0;
    float bar1 = 1000000000.0;
    ivec2 shift[9] = {
            ivec2(0,0),
                ivec2(-1,1), ivec2(0,1), ivec2(1,1),     
                
                ivec2(-1, -1), ivec2(0,-1), ivec2(1,-1),
                ivec2(1,0), ivec2(-1,0)
                };
   

    float r_ = gridsize/2.0-60;
    for (int i=0 ; i<1; i++)
    {
        int x = x0 + shift[i].x;
        int y = y0 + shift[i].y;
        int cellID = (y * int(gridMaxX)) + x;
        if (x<0 || x> gridMaxX) continue;
        if (cellID<indexMax && cellID>=0) {
            float height2_ = popu[cellID].popu;

            float xcen, ycen;
            float r = gridsize/7.5;
            xcen = x*gridsize + float(basexy.x) + gridsize/2.0;
            ycen = y*gridsize + float(basexy.y) + gridsize/2.0;
            int ix = (int(xcen) - BASEX) / 100;
			int iy = (int(ycen) - BASEY) / 100;
			int indexThis = iy * (MAXX - BASEX)/100 + ix;
     

             if (height2_>0) {
                float hNormal =  height2_/4200.0f; //일단 12시로		
	            float h = hNormal * scaleZ * 200; 
                float multipleV = 1.0f;
                float mixC = clamp(hNormal*multipleV,0,1);

                float d0 = sdCappedCylinder(pos - vec3(xcen, ycen, 0), h, r_+30); 
                float d1 = sdCappedCylinder(pos - vec3(xcen, ycen, 0), h+30, r_);
                float d =  opSubtraction(d1, d0);
                bar1 = min(max(d,0.0), bar1);
             }


        }


    }
    
    float judge = 1.0;

    vec2 res;
    float bar;
    if (bar0<bar1) {
        res.y = 1.0 * judge;
        bar = bar0;
    } else {
        res.y = 1.5 * judge;
        bar = bar1;
    }
    
    float ground = pos.z;
    
    res.x = min(ground, bar);
    if (ground<=bar) res.y = 2.0;

    return res;

}




float calcOcclusion( in vec3 pos, in vec3 nor)
{
	float occ = 0.0;
    float sca = 0.0008; //공간 스케일에 따라 다르게 준다.
    for( int i=0; i<5; i++ )
    {
        float h = (worldScale*10.0f) + (worldScale*300.0f)*float(i)/4.0; //공간 스케일에 따라 다르게 준다.
        vec3 opos = pos + h*nor;
        float d = map( opos).x;
        occ += (h-d)*sca;
        sca *= 0.95;
    }
    return clamp( 1.0 - 2.0*occ, 0.001, 1.0 );
}

// http://iquilezles.org/www/articles/rmshadows/rmshadows.htm
float calcSoftshadow(vec3 ro, vec3 rd, vec3 normal)
{ 
    float tmax = 100000.0;  //공간 스케일에 따라 다르게 준다.
	float res = 1.0;
    float t =  0.01;  //0으로 하면 표면 위에서 에러가 발생한다.
    float ph = 1e20; // big, such that y = 0 on the first iteration
    float k = 2.0; //이 숫자를 낮추면 그림자가 넓게 퍼져서 약간 실사 느낌이 난다.
    //ro += normal * 0.1;

    for(int i=0 ; i<2560 && t<tmax; i++ )
    {
		float h = map( ro + rd*t ).x;
        // improved technique

        // use this if you are getting artifact on the first iteration, or unroll the
        // first iteration out of the loop
        float y = h*h/(2.0*ph);
        float d = pow((h*h-y*y), 0.58);
        //float d = sqrt(h*h-y*y);
        res = min( res, k*d/max(0.0,t-y) );
        ph = h;
       
        if (abs(h)<0.005*t) break; 
        //0으로 하면 안된다. h의 비교기준이 너무 작아지면 거기까지 h가 작아지지 않을 경우에 그냥 밝게 나옴

      //  if( res<0.0001 || t>tmax ) break;   
        t += h * marchingUnit ;
          
    }
    return clamp( res, 0.01, 1.0 );
}

// http://iquilezles.org/www/articles/rmshadows/rmshadows.htm
float calcShadow(vec3 ro, vec3 rd, vec3 normal)
{ 
    float tmax = 100000.0;  //공간 스케일에 따라 다르게 준다.
	float res = 1.0;
    float t =  0.01;  //0으로 하면 표면 위에서 에러가 발생한다.
    float ph = 1e20; // big, such that y = 0 on the first iteration
    float k = 1.3;

    for(int i=0 ; i<1280 && t<tmax; i++ )
    {
		float h = map( ro + rd*t ).x;      
       
        if (abs(h)<0.005*t || t>tmax) return 0.1; 
        //0으로 하면 안된다. h의 비교기준이 너무 작아지면 거기까지 h가 작아지지 않을 경우에 그냥 밝게 나옴
        t += h *0.1;
          
    }
    return clamp( t, 0.01, 1.0 );
}

vec2 castRay( vec3 ro,vec3 rd)
{
    vec2 res = vec2(-1.0,-1.0);   
   
    float tmax = 1000000.5;
	#if 1
    // raytrace bounding plane
    //float tp = (3.5-ro.z)/rd.z;
    //if( tp>0.0 ) tmax = min( tmax, tp );
	#endif    
    
    // raymarch scene
    //float t = hash31(ro.zxy + rd.yzx)*.25;
    float t = 0.1;
    int loop = marchingUnit < 0.05? 15120 : 2560; //자세하면 많이 돌린다.
    for( int i=0; i<loop && t<tmax; i++ )
    {
        vec2 h = map( ro+rd*t);

        if( abs(h.x)<(0.00010*t) )
        {       
            res = vec2(t, h.y); 
            break;
        }
        if (res.y>1.6) t += h.x;// * (pow(0.99, i)); //h.x는 중간 레이에서의 거리
        else t += h.x * marchingUnit;
        //막대들이 겹칠 경우 h.x 보다 작게 더해줘야 하는데, 
        //너무 작게 더하면 i의 루프 수를 올려야 하는 문제가 생긴다.
    }
    
    return res;
}

vec3 calcNormal( vec3 pos)
{
/*
    vec2 e = vec2(1.0,-1.0)*0.5773*(worldScale);
    return normalize( e.xyy*map( pos + e.xyy).x + 
					  e.yyx*map( pos + e.yyx ).x + 
					  e.yxy*map( pos + e.yxy).x + 
					  e.xxx*map( pos + e.xxx).x );
                    
   */ 
    // inspired by tdhooper and klems - a way to prevent the compiler from inlining map() 4 times
    vec3 n = vec3(0.0);
    for( int i=0; i<4; i++ )
    {
        vec3 e = 1000.5773*(2.0*vec3((((i+3)>>1)&1),((i>>1)&1),(i&1))-1.0);
        n += e*map(pos+0.001*e).x;
    }
    return normalize(n);   
   
}




float tri(in float x){return abs(fract(x)-.5);}
vec3 tri3(in vec3 p){return vec3( tri(p.y+tri(p.z*1.)),  tri(p.z+tri(p.x*1.)), tri(p.y+tri(p.x*1.)));}
                                
mat2 m2 = mat2(0.970,  0.242, -0.242,  0.970);

float triNoise3d(in vec3 p, in float spd)
{
    float z=1.4;
	float rz = 0.;
    vec3 bp = p;
	for (float i=0.; i<=3.; i++ )
	{
        vec3 dg = tri3(bp*2.);
        p += (dg+time*0.1*spd);

        bp *= 1.8;
		z *= 1.5;
		p *= 1.2;
        //p.xz*= m2;
        
        rz+= (tri(-p.y+tri(p.x+tri(p.z))))/z;
        bp += 0.14;
	}
	return rz;
}
#define USEVOLUMNFOG 1
float fogmap(in vec3 p, in float d) //이 값이 0.~1.0 사이에서 나와야 함
{
    p.x += time*0.1*1.5;
    p.y += sin(p.x*.5);
    //return triNoise3d( p*2.2/(d+20.), 0.2)   *   (1.-smoothstep(0.0,1.0,p.z))  ;
    #if USEVOLUMNFOG==1
    return triNoise3d( p*2.2/(d+20.), 0.2) * (1.0 - smoothstep(0.0,1.0, p.z));//15));// * smoothstep(0.0,1.0,p.z/6) ;
    #else
    return triNoise3d( p*2.2/(d+20.), 0.2) * (1.0 - smoothstep(0.0,1.0, p.z/15)) * smoothstep(0.0,1.0,p.z/6) ;
    #endif
}

float scaleFog = 0.001;

vec3 fog(in vec3 col, in vec3 ro, in vec3 rd, in float mt)
{
    vec3 roRe = ro;
    vec3  sun_lig = normalize(sunlight);
    float mul = 1.0;
    ro -= vec3(1074000., 1728000., .0);
    ro *= scaleFog;
    mt *= scaleFog;
    float d = 0.5; //시선으로부터 이만큼 거리 후에 안개 시작  
  #if USEVOLUMNFOG==1
    for(int i=0; i<550; i++)
  #else
    for(int i=0; i<450; i++)
  #endif
    {
        vec3  pos = ro + rd*d;
        
        #if USEVOLUMNFOG==1
        float shadow =1.0;
        {  
            float dRe = d / scaleFog;
            shadow= clamp(calcSoftshadow(roRe + rd*dRe, sun_lig, vec3(1)),0,1);
        }
        #endif
        float rz0 = fogmap(pos, d); //현재 위치에서 fogmap
        float rz1 = fogmap(pos+ (0.8 - float(i)*0.1) ,  d);

		float grd =  clamp((rz0 - rz1 )*3., 0.1, 1. ); //이 값이 0.1~1 사이에서 나오게 됨

        vec3 col2 = (   vec3(.9,0.7,.3)*.5 + .5*vec3(.5, .8, 1.) * (1.7-grd )  ) *0.55; // ok.
        //vec3 col2 = vec3(1,0,0);
        //col = mix( col, col2,  clamp( rz0* smoothstep( d-0.4, d + 2.0 + d *.75 ,  (mt-d)*50) ,  0.,  1. ) ); //ok   
       
       
        
        #if USEVOLUMNFOG==1
        col = mix( col, col2,  clamp( rz0* pow(clamp((mt-d),0.0,1.0),3)*2.2*shadow,  0.0,  1.0 ) ); //ok     
        d += 0.2;
        #else
        col = mix( col, col2,  clamp( rz0* pow(clamp((mt-d),0.0,1.0),3)*2.2,  0.0,  1.0 ) ); //ok     
        d *= 1.2;
        #endif
        if (d>mt)break;
    }
    return col;
}


vec3 render(  in vec3 ro, in vec3 rd, out vec3 pos, out vec3 ref)
{ 
    float focc = 1.0;
    float ks = 0.35; //specular 에 쓰임 매트하고 반짝이는 정도
    vec3 col = vec3(0.1);

    vec2 res = castRay(ro.xyz, rd); //256 maps
    float t = res.x;
    pos = ro + t*rd;

    float depthPrevious = texelFetch(sDepth, ivec2(gl_FragCoord.xy),0).r;
    vec4 fragcoord = projection * modelview * vec4(pos.xyz, 1.0);
    vec4 fragcoord0 =projection * modelview * vec4(pos.xy, 0, 1.0);
    float depthNow = fragcoord.z/fragcoord.w *0.5 + 0.5;
    float depthNow0 = fragcoord0.z/fragcoord0.w;
    
    if (true) //여기서 depth test를 해버리면, 바닥의 경우에 기존 렌더링과 겹치는 부분으로 인해 일렁거림이 생긴다.
    //if (depthPrevious >=  depthNow)
    //if (depthNow0 > depthNow) 
    {
        gl_FragDepth = depthNow;
        vec3 nor = calcNormal( pos); //4 maps
        ref = reflect( rd, nor );
        
        float occ = calcOcclusion( pos, nor)*focc; // 5 maps
        float fre = clamp(1.0+dot(nor,rd),0.0,1.0);

        vec3  sun_lig = normalize(sunlight);
        float sun_dif = clamp(dot( nor, sun_lig ), 0.0, 1.0 );
        vec3  sun_hal = normalize( sun_lig-rd );
        float sun_sha = clamp(calcSoftshadow( pos, sun_lig, nor),0.1,1.0); //32 maps
        float sun_spe = ks*pow(clamp(dot(nor,sun_hal),0.0,1.0),8.0)*sun_dif*(0.04+0.96*pow(clamp(1.0+dot(sun_hal,rd),0.0,1.0),5.0));
		
        float sky_dif = sqrt(clamp( 0.5+0.5*nor.z, 0.0, 1.0 ));
        float sky_spe = ks*smoothstep( 0.0, 0.5, ref.z )*(0.04+0.96*pow(fre,4.0));

        //pos.z 는 실제 스케일에 따라서 다른 값을 곱해주어야 한다.
        float bou_dif = sqrt(clamp( 0.1-0.9*nor.z, 0.0, 1.0 ))*clamp(1.0-0.0005*pos.z,0.0,1.0);
        float bac_dif = clamp(0.1+0.9*dot( nor, normalize(vec3(-sun_lig.x,-sun_lig.y,0.0))), 0.0, 1.0 );
        float sss_dif = fre*sky_dif*(0.25+0.75*sun_dif*sun_sha);

        float rim_power = 2.9; //커지면 rim 이 얇아진다. 작으면 두꺼워진다.
        float rim_color =pow( smoothstep(0.0,1.0, 1.0-dot(nor, -rd)), rim_power)* (0.5*(1.0+dot(rd, sunlight)));

        if (res.y > 1.9)
        {
            col = groundRendered.xyz *0.024*  vec3(1.9, 4.0, 5.8);

        } else 
        {
            vec3 objectColor;// = mix(materialColor1, materialColor2, clamp(mixC,0,1));
            if (res.y==1.5) objectColor = materialColor2;
            else if (res.y==1.0 ) objectColor = materialColor1;
            else if (res.y== -1.5) objectColor = materialColor3;
            else if (res.y== -1.0) objectColor = materialColor4;
            //col = 0.007*  (vec3(1.9, 4.0, 5.8) + pow(mixC,0.5) * 22.5* objectColor);
            col = 0.04 *  vec3(1.9, 4.0, 5.8) * objectColor;
            
        }

        

        vec3 lin = vec3(0.0);
        lin += sun_dif*vec3(8.10,6.00,4.20)*vec3(sun_sha, sun_sha*sun_sha*0.5+0.5*sun_sha, sun_sha*sun_sha);
        lin += sky_dif*vec3(0.50,0.70,1.00)*occ;
        lin += bou_dif*vec3(0.20,0.70,0.10)*occ;
        lin += bac_dif*vec3(0.45,0.35,0.25)*occ;
        lin += 1*sss_dif*vec3(3.25,2.75,2.50)*occ;
        lin += 10*rim_color*vec3(0.9, 0.5, 0.1);

        col = col*lin;
        col += sun_spe*vec3(9.90,8.10,6.30)*sun_sha;
        col += sky_spe*vec3(0.20,0.30,0.65)*occ*occ;
        col = pow(col,vec3(0.8,0.9,1.0) );
        
        
        // fog00
        float fogT =t;
       // if (false)
        if (pos.z<1000) //안개 높이 설정
        {
            float halfFog = smoothstep(1000,400,pos.z);
            vec3 result = fog(col, ro, rd, t);
         #if USEVOLUMNFOG==1
            col = mix(col, result, halfFog);//*smoothstep(1,0,abs(sun_lig.z)));
         #else
            col = mix(col, result, halfFog*smoothstep(1,0,abs(sun_lig.z)));
         #endif
            
         }
        
        //col = fog(col, ro, rd, t);

        //col = mix( col, vec3(0.9,0.7,0.5), 1.0-exp( -pow(0.000013*fogT, 3)) );

        return col;
       // return res.y>0? col: vec3(0);
    }
    else {
        gl_FragDepth = depthPrevious;
        return vec3(0);

    }
   
}



void main()
{
    vec3 tot = vec3(0.0);
 
    groundRendered = texelFetch(sColor , ivec2(gl_FragCoord.xy),0);

    

    vec2 p = (-iResolution + 2.0*gl_FragCoord.xy)/iResolution;
  
    vec4 poScreen = pmvInversed * vec4(p, 0, 1.0);   
    poScreen /= poScreen.w;
    vec4 ro_ = pmvInversed * vec4(0, 0, -1, 0);
    ro_ /= ro_.w;

    // ray origin
    vec3 ro = ro_.xyz;
    // ray direction
    vec3 rd = normalize(poScreen.xyz - ro);

    vec3 pos, ref;
  
    vec3 col = render( ro, rd, pos, ref);
  
    //이렇게 하면 1차 반사가 된다.
    //vec3 pos1, ref1;
    //vec3 refcol = render( pos, ref, pos1, ref1);
    //col = mix(col, refcol, 0.5);

    // color grading
    col = col*vec3(1.11,0.89,0.79);

    // compress        
    col = 1.35*col/(1.0+col);
        
    // gamma
    col = pow( col, vec3(0.4545) );

    tot += col;


    // s-surve    
    tot = clamp(tot,0.0,1.0);
    tot = tot*tot*(3.0-1.5*tot);

    // vignetting        
    vec2 q = gl_FragCoord.xy/iResolution.xy;
    tot *= 0.5 + 0.5*pow(16.0*q.x*q.y*(1.0-q.x)*(1.0-q.y),0.25);
  
    // output    
    fragColor = vec4( tot, 1.0 );
  
}