/*
Read
http://z0b.kapsi.fi/snippets.php
before using this code. Thank you.
*/
/*
Creates a "shade map". Takes a sun (or any light, actually) position in
latitude/longitude coordinates and the image size, then creates an 8-bit image
showing light/dark areas and smooth transitions between them.
Latitude (Y) is between -90.0 and +90.0 and longitude (X) is between -180.0 and +180.0.
The 'minShade' is the minimum light level. It can be used to prevent "too dark"
dark areas. It must be between 0 and 1, but values > 0.5 usually don't look good.
The output is an 8-bit transition image. It can be used directly to blend two
images, for example. 255 is fully lit, 0 is fully dark.
Returns false if fails, true if OK.
*/
static bool shadeMap(const float sunLatitude, const float sunLongitude, // sun (light) position
const uint32_t width, const uint32_t height, // output size
std::vector &output, // output
const float minShade = 0.0f) // darkest shade
{
if (width < 1 || height < 1 || minShade < 0.0f || minShade > 1.0f)
return false;
// converts from degrees to radians (always multiply)
static const float kToRad = M_PI / 180.0f;
// This specifies how far away the light source is. If it's too small, the
// transitions become too fuzzy. If it's too large, the transitions become
// binary (dark or light, nothing between). We could use real Earth and Sun
// diameters and distance between the two, then fudge the numbers, but why
// bother? The distance varies anyway and the Sun isn't a point-light source.
// Experiments have shown that values between 5 and 15 work well and
// produces nice and reasonable results. Feel free to adjust.
static const float kSunDistance = 10.0f;
// use spherical mapping to project the 2D point into 3D space
const float sunPos[3] =
{
sinf(sunLongitude * kToRad) * cosf(sunLatitude * kToRad) * kSunDistance,
sinf( sunLatitude * kToRad) * kSunDistance,
cosf(sunLongitude * kToRad) * cosf(sunLatitude * kToRad) * kSunDistance
};
// We'll project the 2D map points into 3D space in the same way as the sun
// vector. Because the X is same for every column on every row, we can
// precalculate sin(x) and cos(x) values. It saves us a LOT of repetitive
// calculations.
float sinX[width],
cosX[width];
for (uint32_t x = 0; x < width; x++) {
// convert X to longitude, then to radians
const float v = ((float)x / (float)width - 0.5f) * 360.0f * kToRad;
sinX[x] = sinf(v);
cosX[x] = cosf(v);
}
output.resize(width * height);
uint8_t *out = &output[0];
for (uint32_t y = 0; y < height; y++) {
// Again, use spherical mapping to project the 2D points,
// this gives us surface normals. We're doing a lot of dot
// products and most of the values don't change, so why not
// precompute them? We've precomputed the sin(x) and cos(x)
// values already.
// Convert Y to latitude, then to radians. Project the constants
// (for this row) into 3D space and precalculate every dot
// product factor we can.
const float lat = ((float)y / (float)height - 0.5f) * -180.0f * kToRad,
latS = sinf(lat) * sunPos[1],
latC0 = cosf(lat) * sunPos[0],
latC2 = cosf(lat) * sunPos[2];
for (uint32_t x = 0; x < width; x++) {
// Dot product the surface normal with the sun vector.
// It might not look like it, but the next line actually
// combines a dot product and sunPos[] -style projection.
// We just have precalculated most of the numbers already.
const float dot = sinX[x] * latC0 + cosX[x] * latC2 + latS;
// negative dot products mean the point is "behind" the
// sphere surface and thus receives no light at all
*out++ = (unsigned)(min(max(dot, minShade), 1.0f) * 255.0f) & 0xFF;
}
}
return true;
}