se.llbit.chunky.renderer.scene.Sun.java Source code

Java tutorial

Introduction

Here is the source code for se.llbit.chunky.renderer.scene.Sun.java

Source

/* Copyright (c) 2012-2014 Jesper qvist <jesper@llbit.se>
 *
 * This file is part of Chunky.
 *
 * Chunky is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Chunky is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * You should have received a copy of the GNU General Public License
 * along with Chunky.  If not, see <http://www.gnu.org/licenses/>.
 */
package se.llbit.chunky.renderer.scene;

import java.util.Random;

import org.apache.commons.math3.util.FastMath;

import se.llbit.chunky.renderer.Refreshable;
import se.llbit.chunky.resources.Texture;
import se.llbit.json.JsonObject;
import se.llbit.math.QuickMath;
import se.llbit.math.Ray;
import se.llbit.math.Vector3d;
import se.llbit.util.JSONifiable;

/**
 * Sun model for ray tracing
 * @author Jesper qvist <jesper@llbit.se>
 */
public class Sun implements JSONifiable {

    /**
     * Default sun intensity
     */
    public static final double DEFAULT_INTENSITY = 1.25;

    /**
     * Maximum sun intensity
     */
    public static final double MAX_INTENSITY = 50;

    /**
     * Minimum sun intensity
     */
    public static final double MIN_INTENSITY = 0.1;

    private static final double xZenithChroma[][] = { { 0.00166, -0.00375, 0.00209, 0 },
            { -0.02903, 0.06377, -0.03203, 0.00394 }, { 0.11693, -0.21196, 0.06052, 0.25886 }, };
    private static final double yZenithChroma[][] = { { 0.00275, -0.00610, 0.00317, 0 },
            { -0.04214, 0.08970, -0.04153, 0.00516 }, { 0.15346, -0.26756, 0.06670, 0.26688 }, };
    private static final double mdx[][] = { { -0.0193, -0.2592 }, { -0.0665, 0.0008 }, { -0.0004, 0.2125 },
            { -0.0641, -0.8989 }, { -0.0033, 0.0452 } };
    private static final double mdy[][] = { { -0.0167, -0.2608 }, { -0.0950, 0.0092 }, { -0.0079, 0.2102 },
            { -0.0441, -1.6537 }, { -0.0109, 0.0529 } };
    private static final double mdY[][] = { { 0.1787, -1.4630 }, { -0.3554, 0.4275 }, { -0.0227, 5.3251 },
            { 0.1206, -2.5771 }, { -0.0670, 0.3703 } };

    private static double turb = 2.5;
    private static double turb2 = turb * turb;
    private static Vector3d A = new Vector3d();
    private static Vector3d B = new Vector3d();
    private static Vector3d C = new Vector3d();
    private static Vector3d D = new Vector3d();
    private static Vector3d E = new Vector3d();

    /**
     * Sun texture
     */
    public static Texture texture = new Texture();

    static {
        A.x = mdx[0][0] * turb + mdx[0][1];
        B.x = mdx[1][0] * turb + mdx[1][1];
        C.x = mdx[2][0] * turb + mdx[2][1];
        D.x = mdx[3][0] * turb + mdx[3][1];
        E.x = mdx[4][0] * turb + mdx[4][1];

        A.y = mdy[0][0] * turb + mdy[0][1];
        B.y = mdy[1][0] * turb + mdy[1][1];
        C.y = mdy[2][0] * turb + mdy[2][1];
        D.y = mdy[3][0] * turb + mdy[3][1];
        E.y = mdy[4][0] * turb + mdy[4][1];

        A.z = mdY[0][0] * turb + mdY[0][1];
        B.z = mdY[1][0] * turb + mdY[1][1];
        C.z = mdY[2][0] * turb + mdY[2][1];
        D.z = mdY[3][0] * turb + mdY[3][1];
        E.z = mdY[4][0] * turb + mdY[4][1];
    }

    private double zenith_Y;
    private double zenith_x;
    private double zenith_y;
    private double f0_Y;
    private double f0_x;
    private double f0_y;

    private final Refreshable scene;

    /**
     * Sun radius
     */
    public static final double RADIUS = .03;
    @SuppressWarnings("javadoc")
    public static final double RADIUS_COS = FastMath.cos(RADIUS);
    @SuppressWarnings("javadoc")
    public static final double RADIUS_COS_2 = FastMath.cos(RADIUS * 2);
    @SuppressWarnings("javadoc")
    public static final double RADIUS_SIN = FastMath.sin(RADIUS);
    @SuppressWarnings("javadoc")
    public static final double RADIUS_COS_SQ = RADIUS_COS * RADIUS_COS;
    @SuppressWarnings("javadoc")
    public static final double SUN_WEIGHT = 1 - FastMath.sqrt(1 - RADIUS_SIN * RADIUS_SIN);

    private static final double AMBIENT = .3;

    private double intensity = DEFAULT_INTENSITY;

    private double azimuth = Math.PI / 2.5;
    private double altitude = Math.PI / 3;

    // support vectors
    private final Vector3d su = new Vector3d();
    private final Vector3d sv = new Vector3d();

    /**
     * Location of the sun in the sky.
     */
    private final Vector3d sw = new Vector3d();

    protected final Vector3d emittance = new Vector3d(1, 1, 1);

    // final to ensure that we don't do a lot of redundant re-allocation
    private final Vector3d color = new Vector3d(1, 1, 1);

    /**
     * Calculate skylight for ray
     * @param ray
     */
    public void calcSkyLight(Ray ray, double horizonOffset) {
        double cosTheta = ray.d.y;
        cosTheta += horizonOffset * (1 - cosTheta);
        if (cosTheta < 0)
            cosTheta = 0;
        double cosGamma = ray.d.dot(sw);
        double gamma = FastMath.acos(cosGamma);
        double cos2Gamma = cosGamma * cosGamma;
        double x = zenith_x * perezF(cosTheta, gamma, cos2Gamma, A.x, B.x, C.x, D.x, E.x) * f0_x;
        double y = zenith_y * perezF(cosTheta, gamma, cos2Gamma, A.y, B.y, C.y, D.y, E.y) * f0_y;
        double z = zenith_Y * perezF(cosTheta, gamma, cos2Gamma, A.z, B.z, C.z, D.z, E.z) * f0_Y;
        if (y <= Ray.EPSILON) {
            ray.color.set(0, 0, 0, 1);
        } else {
            double f = (z / y);
            double x2 = x * f;
            double y2 = z;
            double z2 = (1 - x - y) * f;
            // Old CIE-to-RGB matrix
            /*ray.color.set(
                  3.2410*x2 - 1.5374*y2 - 0.4986*z2,
                  -0.9692*x2 + 1.8760*y2 + 0.0416*z2,
                  0.0556*x2 - 0.2040*y2 + 1.0570*z2,
                  1);*/
            // new CIE to RGB M^-1 matrix from http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html
            ray.color.set(2.3706743 * x2 - 0.9000405 * y2 - 0.4706338 * z2,
                    -0.513885 * x2 + 1.4253036 * y2 + 0.0885814 * z2,
                    0.0052982 * x2 - 0.0146949 * y2 + 1.0093968 * z2, 1);
            ray.color.scale(0.045);
        }
    }

    private double chroma(double turb, double turb2, double sunTheta, double[][] matrix) {

        double t1 = sunTheta;
        double t2 = t1 * t1;
        double t3 = t1 * t2;

        return turb2 * (matrix[0][0] * t3 + matrix[0][1] * t2 + matrix[0][2] * t1 + matrix[0][3])
                + turb * (matrix[1][0] * t3 + matrix[1][1] * t2 + matrix[1][2] * t1 + matrix[1][3])
                + (matrix[2][0] * t3 + matrix[2][1] * t2 + matrix[2][2] * t1 + matrix[2][3]);
    }

    private static double perezF(double cosTheta, double gamma, double cos2Gamma, double A, double B, double C,
            double D, double E) {

        return (1 + A * FastMath.exp(B / cosTheta)) * (1 + C * FastMath.exp(D * gamma) + E * cos2Gamma);
    }

    /**
     * Create new sun model
     * @param sceneDescription
     */
    public Sun(Refreshable sceneDescription) {
        this.scene = sceneDescription;
        initSun();
    }

    /**
     * Set equal to other sun model
     * @param other
     */
    public void set(Sun other) {
        azimuth = other.azimuth;
        altitude = other.altitude;
        color.set(other.color);
        intensity = other.intensity;
        initSun();
    }

    private void initSun() {
        double theta = azimuth;
        double phi = altitude;

        double r = QuickMath.abs(FastMath.cos(phi));

        sw.set(FastMath.cos(theta) * r, FastMath.sin(phi), FastMath.sin(theta) * r);

        if (QuickMath.abs(sw.x) > .1) {
            su.set(0, 1, 0);
        } else {
            su.set(1, 0, 0);
        }
        sv.cross(sw, su);
        sv.normalize();
        su.cross(sv, sw);

        emittance.set(color);
        emittance.scale(FastMath.pow(intensity, Scene.DEFAULT_GAMMA));

        updateSkylightValues();
    }

    /**
     * Angle of the sun around the horizon, measured from north
     * @param value
     */
    public void setAzimuth(double value) {
        azimuth = QuickMath.modulo(value, Math.PI * 2);
        initSun();
        scene.refresh();
    }

    /**
     * Sun altitude from the horizon
     * @param value
     */
    public void setAltitude(double value) {
        altitude = QuickMath.clamp(value, 0, Math.PI / 2);
        initSun();
        scene.refresh();
    }

    /**
     * @return Zenith angle
     */
    public double getAltitude() {
        return altitude;
    }

    /**
     * @return Azimuth
     */
    public double getAzimuth() {
        return azimuth;
    }

    /**
     * Check if the ray intersects the sun
     * @param ray
     * @return <code>true</code> if the ray intersects the sun model
     */
    public boolean intersect(Ray ray) {
        /*double dot = ray.d.x * sw.x + ray.d.y * sw.y + ray.d.z * sw.z;
        if (dot >= RADIUS_COS) {
           ray.color.x = emittance.x * 3;
           ray.color.y = emittance.y * 3;
           ray.color.z = emittance.z * 3;
           ray.hit = true;
           return true;
        }*/

        if (ray.d.dot(sw) < .5)
            return false;

        double WIDTH = RADIUS * 4;
        double WIDTH2 = WIDTH * 2;
        double a;
        a = Math.PI / 2 - FastMath.acos(ray.d.dot(su)) + WIDTH;
        if (a >= 0 && a < WIDTH2) {
            double b = Math.PI / 2 - FastMath.acos(ray.d.dot(sv)) + WIDTH;
            if (b >= 0 && b < WIDTH2) {
                texture.getColor(a / WIDTH2, b / WIDTH2, ray.color);
                ray.color.x *= emittance.x * 10;
                ray.color.y *= emittance.y * 10;
                ray.color.z *= emittance.z * 10;
                return true;
            }
        }

        return false;
    }

    /**
     * Calculate flat shading for ray
     * @param ray
     */
    public void flatShading(Ray ray) {
        double shading = ray.n.x * sw.x + ray.n.y * sw.y + ray.n.z * sw.z;
        shading = QuickMath.max(AMBIENT, shading);
        ray.color.x *= emittance.x * shading;
        ray.color.y *= emittance.y * shading;
        ray.color.z *= emittance.z * shading;
    }

    /**
     * @param newColor
     */
    public void setColor(Vector3d newColor) {
        this.color.set(newColor);
        initSun();
        scene.refresh();
    }

    private void updateSkylightValues() {
        double sunTheta = Math.PI / 2 - altitude;
        double cosTheta = FastMath.cos(sunTheta);
        double cos2Theta = cosTheta * cosTheta;
        double chi = (4.0 / 9.0 - turb / 120.0) * (Math.PI - 2 * sunTheta);
        zenith_Y = (4.0453 * turb - 4.9710) * Math.tan(chi) - 0.2155 * turb + 2.4192;
        zenith_Y = (zenith_Y < 0) ? -zenith_Y : zenith_Y;
        zenith_x = chroma(turb, turb2, sunTheta, xZenithChroma);
        zenith_y = chroma(turb, turb2, sunTheta, yZenithChroma);
        f0_x = 1 / perezF(1, sunTheta, cos2Theta, A.x, B.x, C.x, D.x, E.x);
        f0_y = 1 / perezF(1, sunTheta, cos2Theta, A.y, B.y, C.y, D.y, E.y);
        f0_Y = 1 / perezF(1, sunTheta, cos2Theta, A.z, B.z, C.z, D.z, E.z);
    }

    /**
     * Set the sun intensity
     * @param value
     */
    public void setIntensity(double value) {
        intensity = value;
        initSun();
        scene.refresh();
    }

    /**
     * @return The sun intensity
     */
    public double getIntensity() {
        return intensity;
    }

    /**
     * Point ray in random direction within sun solid angle
     * @param reflected
     * @param random
     */
    public void getRandomSunDirection(Ray reflected, Random random) {
        double x1 = random.nextDouble();
        double x2 = random.nextDouble();
        double cos_a = 1 - x1 + x1 * RADIUS_COS;
        double sin_a = FastMath.sqrt(1 - cos_a * cos_a);
        double phi = 2 * Math.PI * x2;

        Vector3d u = new Vector3d(su);
        Vector3d v = new Vector3d(sv);
        Vector3d w = new Vector3d(sw);

        u.scale(FastMath.cos(phi) * sin_a);
        v.scale(FastMath.sin(phi) * sin_a);
        w.scale(cos_a);

        reflected.d.add(u, v);
        reflected.d.add(w);
        reflected.d.normalize();
    }

    /**
     * Atmospheric extinction and inscattering of ray.
     * @param ray
     * @param s
     * @param attenuation
     */
    public void doAtmos(Ray ray, double s, double attenuation) {
        double Br = 0.00002 * 10;
        double Bm = 0.0007 * 10;
        double g = -.001 * 10000;
        double Fex = FastMath.exp(-(Br + Bm) * s);
        if (attenuation < Ray.EPSILON) {
            ray.color.x *= Fex;
            ray.color.y *= Fex;
            ray.color.z *= Fex;
        } else {
            double theta = ray.d.dot(sw);
            double cos_theta = FastMath.cos(theta);
            double cos2_theta = cos_theta * cos_theta;
            double Brt = (3 / (16 * Math.PI)) * Br * (1 + cos2_theta);
            double Bmt = (1 / (4 * Math.PI)) * Bm * ((1 - g) * (1 - g))
                    / FastMath.pow(1 + g * g + 2 * g * cos_theta, 3 / 2.);
            double Fin = ((Brt + Bmt) / (Br + Bm)) * (1 - Fex);
            ray.color.x = ray.color.x * Fex + attenuation * Fin * emittance.x;
            ray.color.y = ray.color.y * Fex + attenuation * Fin * emittance.y;
            ray.color.z = ray.color.z * Fex + attenuation * Fin * emittance.z;
        }
    }

    /**
     * @param d
     * @return Cosine of angle between sun and vector
     */
    public double theta(Vector3d d) {
        return d.dot(sw);
    }

    private static final double Br = 0.0002;
    private static final double Bm = 0.0009;
    private static final double g = -.0007;

    /**
     * @param s
     * @return Extinction factor
     */
    public double extinction(double s) {
        return FastMath.exp(-(Br + Bm) * s);
    }

    /**
     * @param Fex
     * @param theta
     * @return Inscatter factor
     */
    public double inscatter(double Fex, double theta) {
        double cos_theta = FastMath.cos(theta);
        double cos2_theta = cos_theta * cos_theta;
        double Brt = (3 / (16 * Math.PI)) * Br * (1 + cos2_theta);
        double Bmt = (1 / (4 * Math.PI)) * Bm * ((1 - g) * (1 - g))
                / FastMath.pow(1 + g * g + 2 * g * cos_theta, 3 / 2.);
        return ((Brt + Bmt) / (Br + Bm)) * (1 - Fex);
    }

    @Override
    public JsonObject toJson() {
        JsonObject sun = new JsonObject();
        sun.add("altitude", altitude);
        sun.add("azimuth", azimuth);
        sun.add("intensity", intensity);
        JsonObject colorObj = new JsonObject();
        colorObj.add("red", color.x);
        colorObj.add("green", color.y);
        colorObj.add("blue", color.z);
        sun.add("color", colorObj);
        return sun;
    }

    @Override
    public void fromJson(JsonObject obj) {
        azimuth = obj.get("azimuth").doubleValue(Math.PI / 2.5);
        altitude = obj.get("altitude").doubleValue(Math.PI / 3);
        intensity = obj.get("intensity").doubleValue(DEFAULT_INTENSITY);

        JsonObject colorObj = obj.get("color").object();
        color.x = colorObj.get("red").doubleValue(1);
        color.y = colorObj.get("green").doubleValue(1);
        color.z = colorObj.get("blue").doubleValue(1);
    }

    /**
     * @return sun color
     */
    public Vector3d getColor() {
        return new Vector3d(color);
    }
}