NishitaSky - 图1

    1. // Nishita Sky Adapted for Redshift
    2. // Base Nishita implementation built by Ben Simonds - https://github.com/BenSimonds/NishitaSky/blob/master/NishitaAtmos.osl
    3. /** 
    4. Modified: 2021-03-23 by Nischal Pokhrel for Redshift 3D.
    5. This file is licensed under Apache 2.0 license
    6. Added: Z-Up | Light rotation inputs | Sun, Sky & Mie tint controls | Sun disc + intensity, size & tint control | Exposure, saturation & gamma control
    7. Aces mode | Faux ozone | Stars based on elevation and height.
    8. // Starfield based on shader by Thomas Dinges from https://github.com/sambler/osl-shaders
    9. **/
    11. #include <stdosl.h>
    13. int has_solutions (vector StartPoint, vector UnitVector, float Radius) {
    14.     float b = 2 * dot (UnitVector, StartPoint);
    15.     float bsquared = pow(b,2);
    16.     float c = dot (StartPoint,StartPoint) - pow(Radius,2);
    17.     float fourac = 4 * c;
    19.     if (bsquared > fourac) {
    20.         return 1;
    21.     } else {
    22.         return 0;
    23.     }
    24. }
    26. int linetype (vector StartPoint, vector UnitVector, float Radius) {
    27.     float b = 2 * dot (UnitVector, StartPoint);
    28.     float bsquared = pow(b,2);
    29.     float c = dot (StartPoint,StartPoint) - pow(Radius,2);
    30.     float fourac = 4 * c;
    32.     if (bsquared > fourac) {
    33.         float d1 = (-b + pow(bsquared - fourac, 0.5)) / 2;
    34.         float d2 = (-b - pow(bsquared - fourac, 0.5)) / 2;
    35.         if (d1 < 0 && d2 < 0) {
    36.             return 2;
    37.         } else if (d1 < 0 || d2 < 0) {
    38.             return 1;
    39.         } else {
    40.             return 0;
    41.         }
    42.     } else {
    43.         return 3;
    44.     }
    45. }
    47. vector lmin (vector StartPoint, vector UnitVector, float Radius) {
    48.     float b = 2 * dot (UnitVector, StartPoint);
    49.     float bsquared = pow(b,2);
    50.     float c = dot (StartPoint,StartPoint) - pow(Radius,2);
    51.     float fourac = 4 * c;
    53.     if (bsquared > fourac) {
    54.         float d1 = (-b + pow(bsquared - fourac, 0.5)) / 2;
    55.         float d2 = (-b - pow(bsquared - fourac, 0.5)) / 2;
    56.         return StartPoint + UnitVector * min(d1, d2);
    57.     } else {
    58.         return vector(0);
    59.     }
    60. }
    62. vector lmax (vector StartPoint, vector UnitVector, float Radius) {
    63.     float b = 2 * dot (UnitVector, StartPoint);
    64.     float bsquared = pow(b,2);
    65.     float c = dot (StartPoint,StartPoint) - pow(Radius,2);
    66.     float fourac = 4 * c;
    68.     if (bsquared > fourac) {
    69.         float d1 = (-b + pow(bsquared - fourac, 0.5)) / 2;
    70.         float d2 = (-b - pow(bsquared - fourac, 0.5)) / 2;
    71.         return StartPoint + UnitVector * max(d1, d2);
    72.     } else {
    73.         return vector(0);
    74.     }
    75. }
    77. color softLight (color a, color b){
    78.     color ab = 0.0;
    79.     if (length(b) < 0.5){
    80.         ab = 2*(a*b) + (a*a) * (1-(2*b));
    81.     }
    82.     else {
    83.         ab = sqrt(a) * ((2*b) - 1) + ((2*a) * (1-b));
    84.     }
    85.     return (ab);
    86. }
    88. struct linesample {
    89.     int has_solutions;
    90.     int type;
    91.     vector min;
    92.     vector max;
    93. };
    96. shader nishita_sky 
    97. [[  string help = "Nishita Sky Model for Environments",
    98.     string label = "Nishita Sky" ]]
    99. (
    100.     float Intensity = 20            [[string page = "1: Color Settings"]],
    101.     float  Exposure = 0             [[float slidermin = -16, float slidermax = 16, int connectable = 0, string page = "1: Color Settings"]],
    102.     float  Gamma = 1                 [[float slidermin = 0, float slidermax = 2, int connectable = 0, string page = "1: Color Settings"]],
    103.     float Saturation = 1            [[float slidermin = 0, float slidermax = 1.5, int connectable = 0, string page = "1: Color Settings"]],
    104.     float discIntensity = 10        [[float slidermin = 0, float slidermax = 100, string label = "Disc Intensity", string page = "1: Color Settings"]],
    105.     color discColor = color(1,1,1)    [[string label = "Disc Tint", string page = "1: Color Settings"]],
    106.     color sunTint = color(1, 1, 1)    [[string label = "Sun Halo Tint", string page = "1: Color Settings"]],
    107.     color skyTint = color(1, 1, 1)    [[string label = "Sky Tint", string page = "1: Color Settings"]], 
    108.     int convToAces = 0                [[string widget = "checkBox", string label = "ACES-AP1", int connectable = 0, string page = "1: Color Settings"]],    
    109.     float Scattering = 1            [[float slidermin = 0.1, float slidermax = 1.315, float slidercenter = 1, int connectable = 0, string page = "2: Control Settings"]],
    110.     float ozoneDensity = 0            [[float slidermin = -3, float slidermax = 3, float slidercenter = 0, string label = "Ozone Density Multiplier", string page = "2: Control Settings"]],
    111.     float Height = 200                [[string page = "2: Control Settings"]],
    112.     float discSize = 1                [[float slidermin = 0, float slidermax = 100, string label = "Disc Size", string page = "2: Control Settings"]],
    113.     float Azimuth = 0                [[float slidermin = 0, float slidermax = 360, int connectable = 0, string page = "2: Control Settings"]],
    114.     float Elevation = 0                [[float slidermin = -40, float slidermax = 220, int connectable = 0, string page = "2: Control Settings"]],
    115.     int useLight = 0                [[string widget = "checkBox", string label = "Use External Rotations", string page = "2: Control Settings"]],
    116.     vector Light = vector(0, 0, 0)    [[string label = "Rotation Coordinates", string page = "2: Control Settings"]],
    117.     int flipAxis = 0                [[string widget = "checkBox", string label = "Z-Up", int connectable = 0, string page = "2: Control Settings"]],
    118.     int useStars = 1                [[string widget = "checkBox", string label = "Enable Stars",string page = "3: Stars Settings"]],
    119.     float starsIntensity = 1        [[string label = "Stars Intensity", float slidermin = 0, float slidermax = 4 ,string page = "3: Stars Settings"]],
    120.     float starsFrequency = 0.5        [[string label = "Stars Frequency", float slidermin = 0, float slidermax = 1, string page = "3: Stars Settings"]],
    121.     float starsWidth = 0.5            [[string label = "Stars Width", float slidermin = 0, float slidermax = 1, string page = "3: Stars Settings"]],
    122.     output color Sky = 0.0)
    124. {
    125.     color Rayleigh = 0.0;
    126.     color Mie = 0.0;
    127.     color Ozone = 0.0;
    128.     color outSun = 0.0;
    129.     color Stars = 0.0;
    131.     //Conversion matrices
    132.     matrix yUp = 
    133.     {1,0,0,0,
    134.     0,1,0,0,
    135.     0,0,1,0,
    136.     0,0,0,1};
    137.     matrix zUp = 
    138.     {1,0,0,0,
    139.     0,1*flipAxis,1-flipAxis,0,
    140.     0,1-flipAxis,-1*flipAxis,0,
    141.     0,0,0,1};
    142.     matrix xyz_to_rgb =
    143.     {3.24096994,1.53738318,-0.49861076,0,
    144.     -0.96924364,1.8759675,0.04155506,0,
    145.     0.05563008,-0.20397696,1.05697151,0,
    146.     0,0,0,1
    147.     };
    149.     //Variables:
    150.     float Re = 6360e3; // Radius of the earth
    151.     float Ra = 6460e3; // Radius of the atmosphere
    152.     float Hr = 8000; // Rayleigh scale height
    153.     float Hm = 1200; // Mie scale height
    154.     float Ho = 10000; // Ozone scale height
    155.     float g = 0.76 * clamp(Scattering, 0.1, 1.315); // Anisotropy term for Mie scattering.
    156.     float oD = 0.00000022512612292678; // Ozone Coef at wavelength 18
    158.     vector BetaR = vector(5.8e-6,13.0e-6,22.4e-6);
    159.     vector BetaM_max = vector(9.73e-6);
    160.     vector BetaM_min = vector(13.28e-6);
    161.     vector BetaM =  vector(20e-6);
    162.     vector BetaO = vector(3.808e-6);
    164.     //Clean up inputs:
    165.     vector SunDirection = vector(0, 0, 1);
    166.     SunDirection = normalize(SunDirection);
    167.     vector Direction = normalize(I);
    169.     //Sun azimuth, elevation and misc light logic
    170.     if (useLight == 1){
    171.         if (flipAxis == 0){
    172.             SunDirection = transform(yUp, SunDirection);
    173.             SunDirection = rotate(SunDirection, radians(360) - radians(Light[0]), vector(-1,0,0));
    174.             SunDirection = rotate(SunDirection, radians(Light[1]) + radians(180), vector(0,-1,0));
    175.         }
    177.         if (flipAxis == 1){
    178.             SunDirection = rotate(SunDirection, radians(270) - radians(Light[0]), vector(1,0,0));
    179.             SunDirection = rotate(SunDirection, radians(Light[2]), vector(0,0,-1));        
    180.         }
    181.     }
    182.     else {
    183.         if (flipAxis == 0){
    184.             SunDirection = transform(yUp, SunDirection);
    185.             if (Elevation >= -40){
    186.                 SunDirection = rotate(SunDirection, radians(360) - radians(Elevation), vector(1,0,0));
    187.             }
    188.             if (Azimuth >= 0){
    189.                 SunDirection = rotate(SunDirection, radians(Azimuth), vector(0,1,0));
    190.             }
    191.         }
    193.         if (flipAxis == 1){
    194.             if (Elevation >= -40){
    195.                 SunDirection = rotate(SunDirection, radians(270) - radians(Elevation), vector(1,0,0));
    196.             }
    197.             if (Azimuth >= 0){
    198.                 SunDirection = rotate(SunDirection, radians(Azimuth), vector(0,0,-1));
    199.             }
    200.         }
    201.     }
    203.     //Calcualte Rayleigh and mie phases.
    204.     float mu = dot(Direction, SunDirection);
    205.       float phaseR = (3/(16 * M_PI)) * (1 + mu*mu);
    206.       float phaseM = (3/(8 * M_PI)) * ((1 - g*g) * (1 + mu*mu) / ((2 + g*g) * pow(1 + g*g - 2 * g * mu, 1.5)));
    208.     // Sun Disc
    209.     float sunDiameter = discSize * 0.53 * M_PI / 180.0;
    210.     float sunAmount = step(cos(sunDiameter / 2.0), mu);
    212.     //Stuff that will eventually help form my output.
    213.     vector SumR = vector(0);
    214.     vector SumM = vector(0);
    215.     vector SumO = vector(0);
    217.     //Set up start and end points for my integrals.
    218.     vector Pu = vector(0,0,0);
    219.     vector Pv = vector(0,0,0);
    220.     //Get start and end point for View Line.
    221.     if (Height > 0) {
    222.         vector Pc = vector (0, 0, Re + Height);
    223.         if (flipAxis == 0){
    224.             Pc = transform(yUp, Pc);
    225.             Pc = rotate(Pc, radians(-90), vector(1,0,0));
    226.         }
    227.         if (flipAxis == 1){
    228.             Pc = transform(zUp, Pc);
    229.         }
    231.         Pu = Pc;
    233.         //Ozone
    234.         oD = 0;
    235.         if (ozoneDensity != 0){
    236.             if (Height >= 0.0) {
    237.                 oD = (ozoneDensity * 1000 / -15000.0);
    238.                 oD = clamp(oD, -2000, 3000);
    239.             }            
    240.         }
    241.         //Stars stuff
    242.         point PP = I;
    243.         float val;
    244.         float sf = starsFrequency;
    245.         sf /= 300;
    246.         val = smoothstep(1 - starsWidth, 2, noise("perlin", PP / sf));
    247.         vector cosSunGround = dot(normalize(Pu), normalize(SunDirection));
    248.         vector sunAngle = acos(cosSunGround);
    249.         float sunAngleDeg = degrees(sunAngle[0]);
    250.         float sunAngleDeg_N = (sunAngleDeg - 87) / (135 - 87);
    252.         if (useStars == 1){
    253.             if (Height <= 30000){
    254.                 Stars = val * pow(4, clamp(2.5 * starsIntensity, 0, 4)) * sunAngleDeg_N;
    255.                 Stars =  -sunAmount + Stars;            
    256.             }
    257.             else {
    258.                 Stars = val * pow(1.8, clamp(2.5 * starsIntensity, 0, 4));
    259.             }
    260.             Stars = clamp(Stars, 0, 4);    
    261.         }
    263.         linesample groundline;
    264.         groundline.has_solutions = has_solutions(Pc, Direction, Re);
    265.         groundline.type = linetype(Pc, Direction, Re);
    266.         groundline.min = lmin(Pc, Direction, Re);
    267.         groundline.max = lmax(Pc, Direction, Re);
    269.         linesample atmosphereline;
    270.         atmosphereline.has_solutions = has_solutions(Pc, Direction, Ra);
    271.         atmosphereline.type = linetype(Pc, Direction, Ra);
    272.         atmosphereline.min = lmin(Pc, Direction, Ra);
    273.         atmosphereline.max = lmax(Pc, Direction, Ra);
    275.         Pv = atmosphereline.max;
    276.         if (atmosphereline.has_solutions && atmosphereline.type == 0) {
    277.             //Line starts outside atmoshere, so Pu becomes Line.min
    278.             Pu = atmosphereline.min;
    279.         }
    280.         if (groundline.has_solutions && groundline.type == 0) {
    281.             //Line hits ground, Pv becomes point at which line hits ground, which will be Line.min
    282.             Pv = groundline.min;
    283.             sunAmount = 0;
    284.             Stars = 0;
    285.         }
    286.     }
    288.     // Create samples along view ray.
    289.     int numSamples = 16;
    290.     int numSamplesL = 8;
    292.     float segmentLength = distance (Pu, Pv) / numSamples;
    293.     float opticalDepthR = 0;
    294.       float opticalDepthM = 0;
    295.     float opticalDepthO = 0;
    297.     for (int i = 0; i < numSamples; i++) {
    298.     vector Px = Pu + segmentLength * Direction * (i + 0.5);
    299.     float sampleHeight = length(Px) - Re;
    301.         /* Get optical depth for light at this sample: */
    302.         float Hr_sample = exp(-sampleHeight/Hr) * segmentLength;
    303.         float Hm_sample = exp(-sampleHeight/Hm) * segmentLength;
    304.         float Ho_sample = exp(-sampleHeight/Ho) * segmentLength;
    306.         opticalDepthR += Hr_sample;
    307.         opticalDepthM += Hm_sample;
    308.         opticalDepthO += Ho_sample;
    310.         //Get light depth to sun at this sample.
    311.         float opticalDepthLR = 0;
    312.         float opticalDepthLM = 0;
    313.         float opticalDepthLO = 0;
    315.         linesample atmosphereline_sample;
    316.         atmosphereline_sample.has_solutions = has_solutions(Px, SunDirection, Ra);
    317.         atmosphereline_sample.type = linetype(Px, SunDirection, Ra);
    318.         atmosphereline_sample.min = lmin(Px, SunDirection, Ra);
    319.         atmosphereline_sample.max = lmax(Px, SunDirection, Ra);
    321.         linesample groundline_sample;
    322.         groundline_sample.has_solutions = has_solutions(Px, SunDirection, Re);
    323.         groundline_sample.type = linetype(Px, SunDirection, Re);
    324.         groundline_sample.min = lmin(Px, SunDirection, Re);
    325.         groundline_sample.max = lmax(Px, SunDirection, Re);
    327.         vector Ps = atmosphereline_sample.max;
    328.         int j = 0;
    329.         if (groundline_sample.has_solutions && groundline_sample.type == 0) {
    330.             //Light is below horizon from this point. No sample needed.
    331.             opticalDepthLR += 0;
    332.             opticalDepthLM += 0;
    333.             opticalDepthLO += 0;
    334.         } else {
    335.             float segmentLengthL = distance(Px, Ps) / numSamplesL;
    337.             for (j = 0; j < numSamplesL; j++) {
    338.                 vector Pl = Px + segmentLengthL * SunDirection * (j + 0.5);
    339.                 float sampleHeightL = length(Pl) - Re;
    340.                 if (sampleHeightL < 0) break;
    341.                 // ignore sampleheights < 0 they're inside the planet...
    342.                 float Hlr_sample = exp(-sampleHeightL/Hr)*segmentLengthL;
    343.                 float Hlm_sample = exp(-sampleHeightL/Hm)*segmentLengthL;
    344.                 float Hlo_sample = exp(-sampleHeightL/Ho)*segmentLengthL;
    346.                 opticalDepthLR += Hlr_sample;
    347.                 opticalDepthLM += Hlm_sample;
    348.                 opticalDepthLO += Hlo_sample;
    349.             }
    350.         }
    352.         // With our light samples done we can calculate the attenuation.
    353.         // Only include a sample if we reached j (so not if we hit the break clause for rays that hit the ground in finding the sun)
    354.         if (j == numSamplesL) {
    355.             vector tauR = BetaR * (opticalDepthR + opticalDepthLR);
    356.             vector tauM = BetaM * 1.1 * (opticalDepthM + opticalDepthLM);
    357.             vector tauO = BetaO * (opticalDepthO + opticalDepthLO);
    358.             vector tau = tauR + tauM + tauO;
    359.             vector attnR = vector (exp(-tauR[0]), exp(-tauR[1]), exp(-tauR[2]));
    360.             vector attnM = vector (exp(-tauM));
    361.             vector attenuation = vector (exp(-tau[0]), exp(-tau[1]), exp(-tau[2]));
    363.             SumR += Hr_sample * attenuation;
    364.             SumM += Hm_sample * attenuation; 
    365.             SumO += Ho_sample * attenuation;
    366.         }
    367.     }
    368.     matrix srgbToAcesAP1 = {
    369.      0.6131, 0.0701, 0.0206, 0,
    370.     0.3395, 0.9164, 0.1096, 0,
    371.     0.0474, 0.0135, 0.8698, 0,
    372.     0, 0, 0, 1 };
    373.     Rayleigh = color (SumR * phaseR * BetaR) * Intensity * clamp(skyTint, 0, 1);
    374.     Mie = color (SumM * phaseM * BetaM) * Intensity * clamp(sunTint, 0, 1);
    375.     Ozone = color (SumO * BetaO * oD) * Intensity * clamp(sunTint, 0, 1);
    376.     outSun = color (SumM * phaseM * BetaM * discColor) * sunAmount * pow(2, discIntensity);
    377.     // Subtly merge the rayleigh colors with the stars, cull very dim stars
    378.     // Merge the rest of the atmosphere
    379.     color atmosMix = (softLight((Stars * Stars), Rayleigh) + Rayleigh + Mie + Ozone);
    380.     Sky = (atmosMix + outSun);
    381.     color toHSV = transformc("hsv", Sky);
    382.     toHSV[1] = toHSV[1] * Saturation;
    383.     color toRGB = transformc("hsv", "rgb", toHSV);
    384.     Sky = toRGB;
    385.     Sky = pow(Sky, clamp(Gamma, 0, 2)) * pow(2.0, Exposure);
    386.     float sR = Sky[0];
    387.     float sG = Sky[1];
    388.     float sB = Sky[2];
    389.     if (convToAces == 1){
    390.         Sky = transform(srgbToAcesAP1, vector(sR, sG, sB));
    391.         Sky = pow(Sky, Gamma) * pow(2.0, Exposure);
    392.     }
    393. }