CloudCollectionFS.glsl 10.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263
  1. uniform sampler2D u_noiseTexture;
  2. uniform vec3 u_noiseTextureDimensions;
  3. uniform float u_noiseDetail;
  4. varying vec2 v_offset;
  5. varying vec3 v_maximumSize;
  6. varying vec4 v_color;
  7. varying float v_slice;
  8. varying float v_brightness;
  9. float wrap(float value, float rangeLength) {
  10. if(value < 0.0) {
  11. float absValue = abs(value);
  12. float modValue = mod(absValue, rangeLength);
  13. return mod(rangeLength - modValue, rangeLength);
  14. }
  15. return mod(value, rangeLength);
  16. }
  17. vec3 wrapVec(vec3 value, float rangeLength) {
  18. return vec3(wrap(value.x, rangeLength),
  19. wrap(value.y, rangeLength),
  20. wrap(value.z, rangeLength));
  21. }
  22. float textureSliceWidth = u_noiseTextureDimensions.x;
  23. float noiseTextureRows = u_noiseTextureDimensions.y;
  24. float inverseNoiseTextureRows = u_noiseTextureDimensions.z;
  25. float textureSliceWidthSquared = textureSliceWidth * textureSliceWidth;
  26. vec2 inverseNoiseTextureDimensions = vec2(noiseTextureRows / textureSliceWidthSquared,
  27. inverseNoiseTextureRows / textureSliceWidth);
  28. vec2 voxelToUV(vec3 voxelIndex) {
  29. vec3 wrappedIndex = wrapVec(voxelIndex, textureSliceWidth);
  30. float column = mod(wrappedIndex.z, textureSliceWidth * inverseNoiseTextureRows);
  31. float row = floor(wrappedIndex.z / textureSliceWidth * noiseTextureRows);
  32. float xPixelCoord = wrappedIndex.x + column * textureSliceWidth;
  33. float yPixelCoord = wrappedIndex.y + row * textureSliceWidth;
  34. return vec2(xPixelCoord, yPixelCoord) * inverseNoiseTextureDimensions;
  35. }
  36. // Interpolate a voxel with its neighbor (along the positive X-axis)
  37. vec4 lerpSamplesX(vec3 voxelIndex, float x) {
  38. vec2 uv0 = voxelToUV(voxelIndex);
  39. vec2 uv1 = voxelToUV(voxelIndex + vec3(1.0, 0.0, 0.0));
  40. vec4 sample0 = texture2D(u_noiseTexture, uv0);
  41. vec4 sample1 = texture2D(u_noiseTexture, uv1);
  42. return mix(sample0, sample1, x);
  43. }
  44. vec4 sampleNoiseTexture(vec3 position) {
  45. vec3 recenteredPos = position + vec3(textureSliceWidth / 2.0);
  46. vec3 lerpValue = fract(recenteredPos);
  47. vec3 voxelIndex = floor(recenteredPos);
  48. vec4 xLerp00 = lerpSamplesX(voxelIndex, lerpValue.x);
  49. vec4 xLerp01 = lerpSamplesX(voxelIndex + vec3(0.0, 0.0, 1.0), lerpValue.x);
  50. vec4 xLerp10 = lerpSamplesX(voxelIndex + vec3(0.0, 1.0, 0.0), lerpValue.x);
  51. vec4 xLerp11 = lerpSamplesX(voxelIndex + vec3(0.0, 1.0, 1.0), lerpValue.x);
  52. vec4 yLerp0 = mix(xLerp00, xLerp10, lerpValue.y);
  53. vec4 yLerp1 = mix(xLerp01, xLerp11, lerpValue.y);
  54. return mix(yLerp0, yLerp1, lerpValue.z);
  55. }
  56. // Intersection with a unit sphere with radius 0.5 at center (0, 0, 0).
  57. bool intersectSphere(vec3 origin, vec3 dir, float slice,
  58. out vec3 point, out vec3 normal) {
  59. float A = dot(dir, dir);
  60. float B = dot(origin, dir);
  61. float C = dot(origin, origin) - 0.25;
  62. float discriminant = (B * B) - (A * C);
  63. if(discriminant < 0.0) {
  64. return false;
  65. }
  66. float root = sqrt(discriminant);
  67. float t = (-B - root) / A;
  68. if(t < 0.0) {
  69. t = (-B + root) / A;
  70. }
  71. point = origin + t * dir;
  72. if(slice >= 0.0) {
  73. point.z = (slice / 2.0) - 0.5;
  74. if(length(point) > 0.5) {
  75. return false;
  76. }
  77. }
  78. normal = normalize(point);
  79. point -= czm_epsilon2 * normal;
  80. return true;
  81. }
  82. // Transforms the ray origin and direction into unit sphere space,
  83. // then transforms the result back into the ellipsoid's space.
  84. bool intersectEllipsoid(vec3 origin, vec3 dir, vec3 center, vec3 scale, float slice,
  85. out vec3 point, out vec3 normal) {
  86. if(scale.x <= 0.01 || scale.y < 0.01 || scale.z < 0.01) {
  87. return false;
  88. }
  89. vec3 o = (origin - center) / scale;
  90. vec3 d = dir / scale;
  91. vec3 p, n;
  92. bool intersected = intersectSphere(o, d, slice, p, n);
  93. if(intersected) {
  94. point = (p * scale) + center;
  95. normal = n;
  96. }
  97. return intersected;
  98. }
  99. // Assume that if phase shift is being called for octave i,
  100. // the frequency is of i - 1. This saves us from doing extra
  101. // division / multiplication operations.
  102. vec2 phaseShift2D(vec2 p, vec2 freq) {
  103. return (czm_pi / 2.0) * sin(freq.yx * p.yx);
  104. }
  105. vec2 phaseShift3D(vec3 p, vec2 freq) {
  106. return phaseShift2D(p.xy, freq) + czm_pi * vec2(sin(freq.x * p.z));
  107. }
  108. // The cloud texture function derived from Gardner's 1985 paper,
  109. // "Visual Simulation of Clouds."
  110. // https://www.cs.drexel.edu/~david/Classes/Papers/p297-gardner.pdf
  111. const float T0 = 0.6; // contrast of the texture pattern
  112. const float k = 0.1; // computed to produce a maximum value of 1
  113. const float C0 = 0.8; // coefficient
  114. const float FX0 = 0.6; // frequency X
  115. const float FY0 = 0.6; // frequency Y
  116. const int octaves = 5;
  117. float T(vec3 point) {
  118. vec2 sum = vec2(0.0);
  119. float Ci = C0;
  120. vec2 FXY = vec2(FX0, FY0);
  121. vec2 PXY = vec2(0.0);
  122. for(int i = 1; i <= octaves; i++) {
  123. PXY = phaseShift3D(point, FXY);
  124. Ci *= 0.707;
  125. FXY *= 2.0;
  126. vec2 sinTerm = sin(FXY * point.xy + PXY);
  127. sum += Ci * sinTerm + vec2(T0);
  128. }
  129. return k * sum.x * sum.y;
  130. }
  131. const float a = 0.5; // fraction of surface reflection due to ambient or scattered light,
  132. const float t = 0.4; // fraction of texture shading
  133. const float s = 0.25; // fraction of specular reflection
  134. float I(float Id, float Is, float It) {
  135. return (1.0 - a) * ((1.0 - t) * ((1.0 - s) * Id + s * Is) + t * It) + a;
  136. }
  137. const vec3 lightDir = normalize(vec3(0.2, -1.0, 0.7));
  138. vec4 drawCloud(vec3 rayOrigin, vec3 rayDir, vec3 cloudCenter, vec3 cloudScale, float cloudSlice,
  139. float brightness) {
  140. vec3 cloudPoint, cloudNormal;
  141. if(!intersectEllipsoid(rayOrigin, rayDir, cloudCenter, cloudScale, cloudSlice,
  142. cloudPoint, cloudNormal)) {
  143. return vec4(0.0);
  144. }
  145. float Id = clamp(dot(cloudNormal, -lightDir), 0.0, 1.0); // diffuse reflection
  146. float Is = max(pow(dot(-lightDir, -rayDir), 2.0), 0.0); // specular reflection
  147. float It = T(cloudPoint); // texture function
  148. float intensity = I(Id, Is, It);
  149. vec3 color = vec3(intensity * clamp(brightness, 0.1, 1.0));
  150. vec4 noise = sampleNoiseTexture(u_noiseDetail * cloudPoint);
  151. float W = noise.x;
  152. float W2 = noise.y;
  153. float W3 = noise.z;
  154. // The dot product between the cloud's normal and the ray's direction is greatest
  155. // in the center of the ellipsoid's surface. It decreases towards the edge.
  156. // Thus, it is used to blur the areas leading to the edges of the ellipsoid,
  157. // so that no harsh lines appear.
  158. // The first (and biggest) layer of worley noise is then subtracted from this.
  159. // The final result is scaled up so that the base cloud is not too translucent.
  160. float ndDot = clamp(dot(cloudNormal, -rayDir), 0.0, 1.0);
  161. float TR = pow(ndDot, 3.0) - W; // translucency
  162. TR *= 1.3;
  163. // Subtracting the second and third layers of worley noise is more complicated.
  164. // If these layers of noise were simply subtracted from the current translucency,
  165. // the shape derived from the first layer of noise would be completely deleted.
  166. // The erosion of this noise should thus be constricted to the edges of the cloud.
  167. // However, because the edges of the ellipsoid were already blurred away, mapping
  168. // the noise to (1.0 - ndDot) will have no impact on most of the cloud's appearance.
  169. // The value of (0.5 - ndDot) provides the best compromise.
  170. float minusDot = 0.5 - ndDot;
  171. // Even with the previous calculation, subtracting the second layer of wnoise
  172. // erode too much of the cloud. The addition of it, however, will detailed
  173. // volume to the cloud. As long as the noise is only added and not subtracted,
  174. // the results are aesthetically pleasing.
  175. // The minusDot product is mapped in a way that it is larger at the edges of
  176. // the ellipsoid, so a subtraction and min operation are used instead of
  177. // an addition and max one.
  178. TR -= min(minusDot * W2, 0.0);
  179. // The third level of worley noise is subtracted from the result, with some
  180. // modifications. First, a scalar is added to minusDot so that the noise
  181. // starts affecting the shape farther away from the center of the ellipsoid's
  182. // surface. Then, it is scaled down so its impact is not too intense.
  183. TR -= 0.8 * (minusDot + 0.25) * W3;
  184. // The texture function's shading does not correlate with the shape of the cloud
  185. // produced by the layers of noise, so an extra shading scalar is calculated.
  186. // The darkest areas of the cloud are assigned to be where the noise erodes
  187. // the cloud the most. This is then interpolated based on the translucency
  188. // and the diffuse shading term of that point in the cloud.
  189. float shading = mix(1.0 - 0.8 * W * W, 1.0, Id * TR);
  190. // To avoid values that are too dark, this scalar is increased by a small amount
  191. // and clamped so it never goes to zero.
  192. shading = clamp(shading + 0.2, 0.3, 1.0);
  193. // Finally, the contrast of the cloud's color is increased.
  194. vec3 finalColor = mix(vec3(0.5), shading * color, 1.15);
  195. return vec4(finalColor, clamp(TR, 0.0, 1.0)) * v_color;
  196. }
  197. void main() {
  198. #ifdef DEBUG_BILLBOARDS
  199. gl_FragColor = vec4(0.0, 0.5, 0.5, 1.0);
  200. #endif
  201. // To avoid calculations with high values,
  202. // we raycast from an arbitrarily smaller space.
  203. vec2 coordinate = v_maximumSize.xy * v_offset;
  204. vec3 ellipsoidScale = 0.82 * v_maximumSize;
  205. vec3 ellipsoidCenter = vec3(0.0);
  206. float zOffset = max(ellipsoidScale.z - 10.0, 0.0);
  207. vec3 eye = vec3(0, 0, -10.0 - zOffset);
  208. vec3 rayDir = normalize(vec3(coordinate, 1.0) - eye);
  209. vec3 rayOrigin = eye;
  210. #ifdef DEBUG_ELLIPSOIDS
  211. vec3 point, normal;
  212. if(intersectEllipsoid(rayOrigin, rayDir, ellipsoidCenter, ellipsoidScale, v_slice,
  213. point, normal)) {
  214. gl_FragColor = v_brightness * v_color;
  215. }
  216. #else
  217. #ifndef DEBUG_BILLBOARDS
  218. vec4 cloud = drawCloud(rayOrigin, rayDir,
  219. ellipsoidCenter, ellipsoidScale, v_slice, v_brightness);
  220. if(cloud.w < 0.01) {
  221. discard;
  222. }
  223. gl_FragColor = cloud;
  224. #endif
  225. #endif
  226. }