package org.codemap.internal;
import static java.lang.Math.PI;
import static java.lang.Math.atan;
import static java.lang.Math.atan2;
import static java.lang.Math.cos;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
import static org.codemap.CodemapCore.colorScheme;
import org.codemap.MapAlgorithm;
import org.codemap.MapSetting;
import org.codemap.util.StopWatch;
/**
* @see http://edndoc.esri.com/arcobjects/9.2/NET/shared/geoprocessing/ spatial_analyst_tools/how_hillshade_works.htm
* @see Burrough, P. A. and McDonell, R.A., 1998. Principles of Geographical Information Systems (Oxford University
* Press, New York), p. 190.
*/
public class ShadeAlgorithm extends MapAlgorithm<double[][]> {
private static final double Z_FACTOR_MULT = 1e-3;
public static final MapSetting<Integer> CONTOUR_STEP = MapSetting.define("CONTOUR_STEP", 10);
@Override
public double[][] call() {
StopWatch stopWatch = new StopWatch("HillShade").start();
// zenith: height of sun over horizon (90 = horizon, 0 = zenith).
double zenithRad = 45 * PI / 180;
// azimuth: direction of sun on x-y-plane
double azimuthRad = (315 - 180) * PI / 180;
// generate magic factor that takes darken factor of color scheme into account
// furthermore, it adapts to the map-size
double z_factor = (1.1 - colorScheme().getDarkenFactor()) * Z_FACTOR_MULT * map.getWidth();
int step = map.get(CONTOUR_STEP);
double[][] hillshade = new double[map.width][map.width];
float[][] DEM = map.getDEM();
int width = map.width;
double contour_darken_factor = colorScheme().getDarkenFactor();
for (int x = 0; x < width; x++) {
for (int y = 0; y < width; y++) {
float here = DEM[x][y];
if (here == 0.0) continue;
// build kernel
int yTop = (y == 0 ? 0 : y - 1);
int xLeft = (x == 0 ? 0 : x - 1);
int yBottom = (y == (width -1) ? (width - 1) : y + 1);
int xRight = (x == (width - 1) ? (width - 1) : x + 1);
float topLeft = DEM[xLeft][yTop];
float top = DEM[x][yTop];
float topRight = DEM[xRight][yTop];
float left = DEM[xLeft][y];
float right = DEM[xRight][y];
float bottomLeft = DEM[xLeft][yBottom];
float bottom = DEM[x][yBottom];
float bottomRight = DEM[xRight][yBottom];
// calculate hillshading
double dx = (topRight + (2 * right) + bottomRight - (topLeft + (2 * left) + bottomLeft)) / 8;
double dy = (bottomLeft + (2 * bottom) + bottomRight - (topLeft + (2 * top) + topRight)) / 8;
double slopeRad = atan(z_factor * sqrt(dx * dx + (dy * dy)));
double aspectRad = atan2(dy, -dx);
double shading = (cos(zenithRad) * cos(slopeRad) + (sin(zenithRad) * sin(slopeRad) * cos(azimuthRad
- aspectRad)));
hillshade[x][y] = shading; //min(0.5, shading);
// calculate contourlines
int topHeight = (int) Math.floor(top / step);
int leftHeight = (int) Math.floor(left / step);
int hereHeight = (int) Math.floor(here / step);
boolean isCoastline = topHeight == 0 || leftHeight == 0 || hereHeight == 0;
boolean isContourLine = (topHeight != hereHeight || leftHeight != hereHeight) && !isCoastline;
if (isContourLine) {
hillshade[x][y] *= contour_darken_factor;
}
}
}
stopWatch.printStop();
return hillshade;
}
}