//    HermiteNoiseマップ作成シェーダ
//        ver 1.0 2020/06/07
//        by yarunashi@dooon


//
//	Code repository for GPU noise development blog
//	http://briansharpe.wordpress.com
//	https://github.com/BrianSharpe
//
//	I'm not one for copyrights.  Use the code however you wish.
//	All I ask is that credit be given back to the blog or myself when appropriate.
//	And also to let me know if you come up with any changes, improvements, thoughts or interesting uses for this stuff. :)
//	Thanks!
//
//	Brian Sharpe
//	brisharpe CIRCLE_A yahoo DOT com
//	http://briansharpe.wordpress.com
//	https://github.com/BrianSharpe
//

void FAST32_hash_3D( 	vec3 gridcell,
                        out vec4 lowz_hash_0,
                        out vec4 lowz_hash_1,
                        out vec4 lowz_hash_2,
                        out vec4 highz_hash_0,
                        out vec4 highz_hash_1,
                        out vec4 highz_hash_2	)		//	generates 3 random numbers for each of the 8 cell corners
{
    //    gridcell is assumed to be an integer coordinate

    //	TODO: 	these constants need tweaked to find the best possible noise.
    //			probably requires some kind of brute force computational searching or something....
    const vec2 OFFSET = vec2( 50.0, 161.0 );
    const float DOMAIN = 69.0;
    const vec3 SOMELARGEFLOATS = vec3( 635.298681, 682.357502, 668.926525 );
    const vec3 ZINC = vec3( 48.500388, 65.294118, 63.934599 );

    //	truncate the domain
    gridcell.xyz = gridcell.xyz - floor(gridcell.xyz * ( 1.0 / DOMAIN )) * DOMAIN;
    vec3 gridcell_inc1 = step( gridcell, vec3( DOMAIN - 1.5 ) ) * ( gridcell + 1.0 );

    //	calculate the noise
    vec4 P = vec4( gridcell.xy, gridcell_inc1.xy ) + OFFSET.xyxy;
    P *= P;
    P = P.xzxz * P.yyww;
    vec3 lowz_mod = vec3( 1.0 / ( SOMELARGEFLOATS.xyz + gridcell.zzz * ZINC.xyz ) );
    vec3 highz_mod = vec3( 1.0 / ( SOMELARGEFLOATS.xyz + gridcell_inc1.zzz * ZINC.xyz ) );
    lowz_hash_0 = fract( P * lowz_mod.xxxx );
    highz_hash_0 = fract( P * highz_mod.xxxx );
    lowz_hash_1 = fract( P * lowz_mod.yyyy );
    highz_hash_1 = fract( P * highz_mod.yyyy );
    lowz_hash_2 = fract( P * lowz_mod.zzzz );
    highz_hash_2 = fract( P * highz_mod.zzzz );
}

//
//	Quintic Hermite Interpolation
//	http://www.rose-hulman.edu/~finn/CCLI/Notes/day09.pdf
//
//  NOTE: maximum value of a hermitequintic interpolation with zero acceleration at the endpoints would be...
//        f(x=0.5) = MAXPOS + MAXVELOCITY * ( ( x - 6x^3 + 8x^4 - 3x^5 ) - ( -4x^3 + 7x^4 -3x^5 ) ) = MAXPOS + MAXVELOCITY * 0.3125
//
//	variable naming conventions:
//	val = value ( position )
//	grad = gradient ( velocity )
//	x = 0.0->1.0 ( time )
//	i = interpolation = a value to be interpolated
//	e = evaluation = a value to be used to calculate the interpolation
//	0 = start
//	1 = end
//
float QuinticHermite( float x, float ival0, float ival1, float egrad0, float egrad1 )		// quintic hermite with start/end acceleration of 0.0
{
    const vec3 C0 = vec3( -15.0, 8.0, 7.0 );
    const vec3 C1 = vec3( 6.0, -3.0, -3.0 );
    const vec3 C2 = vec3( 10.0, -6.0, -4.0 );
    vec3 h123 = ( ( ( C0 + C1 * x ) * x ) + C2 ) * ( x*x*x );
    return ival0 + dot( vec3( (ival1 - ival0), egrad0, egrad1 ), h123.xyz + vec3( 0.0, x, 0.0 ) );
}
vec4 QuinticHermite( float x, vec4 ival0, vec4 ival1, vec4 egrad0, vec4 egrad1 )		// quintic hermite with start/end acceleration of 0.0
{
    const vec3 C0 = vec3( -15.0, 8.0, 7.0 );
    const vec3 C1 = vec3( 6.0, -3.0, -3.0 );
    const vec3 C2 = vec3( 10.0, -6.0, -4.0 );
    vec3 h123 = ( ( ( C0 + C1 * x ) * x ) + C2 ) * ( x*x*x );
    return ival0 + (ival1 - ival0) * h123.xxxx + egrad0 * vec4( h123.y + x ) + egrad1 * h123.zzzz;
}
vec4 QuinticHermite( float x, vec2 igrad0, vec2 igrad1, vec2 egrad0, vec2 egrad1 )		// quintic hermite with start/end position and acceleration of 0.0
{
    const vec3 C0 = vec3( -15.0, 8.0, 7.0 );
    const vec3 C1 = vec3( 6.0, -3.0, -3.0 );
    const vec3 C2 = vec3( 10.0, -6.0, -4.0 );
    vec3 h123 = ( ( ( C0 + C1 * x ) * x ) + C2 ) * ( x*x*x );
    return vec4( egrad1, igrad0 ) * vec4( h123.zz, 1.0, 1.0 ) + vec4( egrad0, h123.xx ) * vec4( vec2( h123.y + x ), (igrad1 - igrad0) );	//	returns vec4( out_ival.xy, out_igrad.xy )
}
void QuinticHermite( 	float x,
                        vec4 ival0, vec4 ival1,			//	values are interpolated using the gradient arguments
                        vec4 igrad_x0, vec4 igrad_x1, 	//	gradients are interpolated using eval gradients of 0.0
                        vec4 igrad_y0, vec4 igrad_y1,
                        vec4 egrad0, vec4 egrad1, 		//	our evaluation gradients
                        out vec4 out_ival, out vec4 out_igrad_x, out vec4 out_igrad_y )	// quintic hermite with start/end acceleration of 0.0
{
    const vec3 C0 = vec3( -15.0, 8.0, 7.0 );
    const vec3 C1 = vec3( 6.0, -3.0, -3.0 );
    const vec3 C2 = vec3( 10.0, -6.0, -4.0 );
    vec3 h123 = ( ( ( C0 + C1 * x ) * x ) + C2 ) * ( x*x*x );
    out_ival = ival0 + (ival1 - ival0) * h123.xxxx + egrad0 * vec4( h123.y + x ) + egrad1 * h123.zzzz;
    out_igrad_x = igrad_x0 + (igrad_x1 - igrad_x0) * h123.xxxx;	//	NOTE: gradients of 0.0
    out_igrad_y = igrad_y0 + (igrad_y1 - igrad_y0) * h123.xxxx;	//	NOTE: gradients of 0.0
}
void QuinticHermite( 	float x,
                        vec4 igrad_x0, vec4 igrad_x1, 	//	gradients are interpolated using eval gradients of 0.0
                        vec4 igrad_y0, vec4 igrad_y1,
                        vec4 egrad0, vec4 egrad1, 		//	our evaluation gradients
                        out vec4 out_ival, out vec4 out_igrad_x, out vec4 out_igrad_y )	// quintic hermite with start/end position and acceleration of 0.0
{
    const vec3 C0 = vec3( -15.0, 8.0, 7.0 );
    const vec3 C1 = vec3( 6.0, -3.0, -3.0 );
    const vec3 C2 = vec3( 10.0, -6.0, -4.0 );
    vec3 h123 = ( ( ( C0 + C1 * x ) * x ) + C2 ) * ( x*x*x );
    out_ival = egrad0 * vec4( h123.y + x ) + egrad1 * h123.zzzz;
    out_igrad_x = igrad_x0 + (igrad_x1 - igrad_x0) * h123.xxxx;	//	NOTE: gradients of 0.0
    out_igrad_y = igrad_y0 + (igrad_y1 - igrad_y0) * h123.xxxx;	//	NOTE: gradients of 0.0
}

//
//	Hermite3D
//	Return value range of -1.0->1.0
//	http://briansharpe.files.wordpress.com/2012/01/hermitesample.jpg
//
float Hermite3D( vec3 P )
{
    //	establish our grid cell and unit position
    vec3 Pi = floor(P);
    vec3 Pf = P - Pi;

    //	calculate the hash.
    //	( various hashing methods listed in order of speed )
    vec4 hash_gradx0, hash_grady0, hash_gradz0, hash_gradx1, hash_grady1, hash_gradz1;
    FAST32_hash_3D( Pi, hash_gradx0, hash_grady0, hash_gradz0, hash_gradx1, hash_grady1, hash_gradz1 );

    //	scale the hash values
    hash_gradx0 = ( hash_gradx0 - 0.49999);
    hash_grady0 = ( hash_grady0 - 0.49999);
    hash_gradz0 = ( hash_gradz0 - 0.49999);
    hash_gradx1 = ( hash_gradx1 - 0.49999);
    hash_grady1 = ( hash_grady1 - 0.49999);
    hash_gradz1 = ( hash_gradz1 - 0.49999);

#if 1
    //	normalize gradients
    vec4 norm0 = inversesqrt( hash_gradx0 * hash_gradx0 + hash_grady0 * hash_grady0 + hash_gradz0 * hash_gradz0 );
    hash_gradx0 *= norm0;
    hash_grady0 *= norm0;
    hash_gradz0 *= norm0;
    vec4 norm1 = inversesqrt( hash_gradx1 * hash_gradx1 + hash_grady1 * hash_grady1 + hash_gradz1 * hash_gradz1 );
    hash_gradx1 *= norm1;
    hash_grady1 *= norm1;
    hash_gradz1 *= norm1;
    const float FINAL_NORM_VAL = 1.8475208614068024464292760976063;
#else
    //	unnormalized gradients
    const float FINAL_NORM_VAL = (1.0/0.46875);  // = 1.0 / ( 0.5 * 0.3125 * 3.0 )
#endif

    //	evaluate the hermite
    vec4 ival_results, igrad_results_x, igrad_results_y;
    QuinticHermite( Pf.z, hash_gradx0, hash_gradx1, hash_grady0, hash_grady1, hash_gradz0, hash_gradz1, ival_results, igrad_results_x, igrad_results_y );
    vec4 qh_results = QuinticHermite( Pf.y, vec4(ival_results.xy, igrad_results_x.xy), vec4(ival_results.zw, igrad_results_x.zw), vec4( igrad_results_y.xy, 0.0, 0.0 ), vec4( igrad_results_y.zw, 0.0, 0.0 ) );
    return QuinticHermite( Pf.x, qh_results.x, qh_results.y, qh_results.z, qh_results.w ) * FINAL_NORM_VAL;
}




// xs_begin
// author : 'yarunashi@dooon'
// arg : { id = '0' name = 'xscale' value = '2' range = '0.0 100.0' step = '0.5'		decimal = '1' }
// arg : { id = '1' name = 'yscale' value = '2' range = '0.0 100.0' step = '0.5'		decimal = '1' }
// arg : { id = '2' name = 'seed' value = '0' range = '0.0 100.0' step = '0.1' decimal = '1' }
// arg : { id = '3' name = 'numcolors' value = '255' range = '1 255' step = '1' decimal = '0' }
// arg : { id = '4' name = 'seamless' value = '0' range = '0 1' step = '1' decimal = '0' }
// arg : { id = '5' name = 'fbm' value = '0' range = '0 1' step = '1' decimal = '0' }
// arg : { id = '6' name = 'hscale' value = '1.0' range = '0.1 100.0' step = '0.1'		decimal = '1' }
// arg : { id = '7' name = 'thickness' value = '1' range = '1 255' step = '1' decimal = '0' }
// xs_end

float xscale = (i_args[0] == 0) ? 2.0 : iArgs[0] ; 
float yscale = (i_args[1] == 0) ? 2.0 : iArgs[1] ; 
float seed = (i_args[2] == 0) ? 0.0 : iArgs[2] ;
float step = (i_args[3] == 0) ? 255.0 : iArgs[3] ;
float seamless = (i_args[4] == 0) ? 0.0 : 1.0;
float fbm = (i_args[5] == 0) ? 0.0 : 1.0;
float hscale = (iArgs[6] == 0) ? 1.0 : iArgs[6] ; 
float thickness = (iArgs[7] == 0) ? 1.0 : iArgs[7] ; 


// シームありノイズ生成
float smnoise(vec2 v)
{
    vec3 p = vec3(v.x * xscale / i_volume_size.x, v.y * yscale / i_volume_size.y, seed);
    float r = Hermite3D(p);
    return (r + 1.0) / 2.0; //0.0~1.0
}

//フラクショナルブラウン運動
#define NUM_OCTAVES 5
float fsmnoise(vec2 x) {
	float v = 0.0;
	float a = 0.5;
	vec2 shift = vec2(100);
	// Rotate to reduce axial bias
    mat2 rot = mat2(cos(0.5), sin(0.5), -sin(0.5), cos(0.50));
	for (int i = 0; i < NUM_OCTAVES; ++i) {
		v += a * smnoise(x);
		x = rot * x * 2.0 + shift;
		a *= 0.5;
	}
	return v;
}


// シームレスノイズ生成
float smlessnoise(vec2 v){
	vec2 p = mod(v, i_volume_size.xy);
    vec2 q = p / i_volume_size.xy;
    vec2 r = i_volume_size.xy;
    return smnoise(vec2(p.x,       p.y      )) *        q.x  *        q.y  +
           smnoise(vec2(p.x,       p.y + r.y)) *        q.x  * (1.0 - q.y) +
           smnoise(vec2(p.x + r.x, p.y      )) * (1.0 - q.x) *        q.y  +
           smnoise(vec2(p.x + r.x, p.y + r.y)) * (1.0 - q.x) * (1.0 - q.y);
}

float fsmlessnoise(vec2 v){
	vec2 p = mod(v, i_volume_size.xy);
    vec2 q = p / i_volume_size.xy;
    vec2 r = i_volume_size.xy;
    return fsmnoise(vec2(p.x,       p.y      )) *        q.x  *        q.y  +
           fsmnoise(vec2(p.x,       p.y + r.y)) *        q.x  * (1.0 - q.y) +
           fsmnoise(vec2(p.x + r.x, p.y      )) * (1.0 - q.x) *        q.y  +
           fsmnoise(vec2(p.x + r.x, p.y + r.y)) * (1.0 - q.x) * (1.0 - q.y);
}

float makemap(vec3 v) {  
    float r = 0.0;
    if (seamless >= 1.0) {
        // seamless noise
        if (fbm >= 1.0) {
            r = fsmlessnoise(v.xy);
        } else {
            r = smlessnoise(v.xy);
        }
    } else {
        if (fbm >= 1.0) {
            r = fsmnoise(v.xy);
        } else {
            r = smnoise(v.xy);
        }
    }
    r = r * 255;
    return r - mod(r, 255.0 / step) + 1;
}

float map(vec3 v) {  

	//マップの色1-255を取得
	float index = makemap(v);

	//色の高さを反転
	float height = 256.0 - index;

	//空間の高さに合わせる
	height = (height / 255.0) * i_volume_size.z * hscale;

    if (i_volume_size.z > 1.0) {
        //色の幅を調整
        index = 256.0 - height;
        index = index - mod(index, thickness) + 1;
    }

    return ( floor(v.z) < height ? index : 0.0 );
}