/* 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;
import se.llbit.util.VectorPool;
/**
* 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;
ray.hit = true;
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
* @param vectorPool
*/
public void getRandomSunDirection(Ray reflected, Random random, VectorPool vectorPool) {
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 = vectorPool.get(su);
Vector3d v = vectorPool.get(sv);
Vector3d w = vectorPool.get(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();
vectorPool.dispose(u);
vectorPool.dispose(v);
vectorPool.dispose(w);
}
/**
* 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);
}
}