1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/terrain.cc	Tue Feb 25 23:47:20 2014 +0200
     1.3 @@ -0,0 +1,349 @@
     1.4 +#include <goat3dgfx/goat3dgfx.h>
     1.5 +#include "terrain.h"
     1.6 +#include "imago2.h"
     1.7 +
     1.8 +using namespace goatgfx;
     1.9 +
    1.10 +Terrain::Terrain()
    1.11 +{
    1.12 +	xsz = ysz = 0;
    1.13 +	dx = dy = 0;
    1.14 +	hmap = 0;
    1.15 +	scale = height_scale = 1;
    1.16 +
    1.17 +	dlist = dlist_norm = 0;
    1.18 +
    1.19 +	dlist_sub[0] = dlist_sub[1] = dlist_norm_sub[0] = dlist_norm_sub[1] = -1;
    1.20 +	dlist_norm_sz = -1.0;
    1.21 +}
    1.22 +
    1.23 +Terrain::~Terrain()
    1.24 +{
    1.25 +	if(hmap) {
    1.26 +		img_free_pixels(hmap);
    1.27 +	}
    1.28 +
    1.29 +	glDeleteLists(dlist, 2);
    1.30 +}
    1.31 +
    1.32 +void Terrain::set_scale(float s)
    1.33 +{
    1.34 +	if(s != scale) {
    1.35 +		scale = s;
    1.36 +		invalidate_mesh();
    1.37 +	}
    1.38 +}
    1.39 +
    1.40 +float Terrain::get_scale() const
    1.41 +{
    1.42 +	return scale;
    1.43 +}
    1.44 +
    1.45 +void Terrain::set_height_scale(float s)
    1.46 +{
    1.47 +	if(s != height_scale) {
    1.48 +		height_scale = s;
    1.49 +		invalidate_mesh();
    1.50 +	}
    1.51 +}
    1.52 +
    1.53 +float Terrain::get_height_scale() const
    1.54 +{
    1.55 +	return height_scale;
    1.56 +}
    1.57 +
    1.58 +bool Terrain::load(const char *fname)
    1.59 +{
    1.60 +	if(!(hmap = (float*)img_load_pixels(fname, &xsz, &ysz, IMG_FMT_GREYF))) {
    1.61 +		fprintf(stderr, "failed to load image: %s\n", fname);
    1.62 +		return false;
    1.63 +	}
    1.64 +
    1.65 +	dx = 1.0 / (float)(xsz - 1);
    1.66 +	dy = 1.0 / (float)(ysz - 1);
    1.67 +	return true;
    1.68 +}
    1.69 +
    1.70 +Vector2 Terrain::world_to_uv(const Vector2 &pt) const
    1.71 +{
    1.72 +	return pt  / scale + 0.5;
    1.73 +}
    1.74 +
    1.75 +Vector2 Terrain::uv_to_world(const Vector2 &pt) const
    1.76 +{
    1.77 +	return (pt - 0.5) * scale;
    1.78 +}
    1.79 +
    1.80 +float Terrain::lookup(int x, int y) const
    1.81 +{
    1.82 +	int px = x < 0 ? 0 : (x >= xsz ? x = xsz - 1 : x);
    1.83 +	int py = y < 0 ? 0 : (y >= ysz ? y = ysz - 1 : y);
    1.84 +	return hmap[py * xsz + px];
    1.85 +}
    1.86 +
    1.87 +float Terrain::get_height(const Vector2 &pt) const
    1.88 +{
    1.89 +	float fxsz = (float)xsz;
    1.90 +	float fysz = (float)ysz;
    1.91 +
    1.92 +	// the floor of x * xsz is the pixel x coordinate (e.g. [0.0, 1.0, 2.0, ... 127.0] for 128x128)
    1.93 +	float floorx = floor(pt.x * fxsz);
    1.94 +	float floory = floor(pt.y * fysz);
    1.95 +
    1.96 +	// (x0, y0) is the lower left pixel in normalized coords of the four pixels surrounding pt
    1.97 +	float x0 = floorx / fxsz;
    1.98 +	float y0 = floory / fysz;
    1.99 +
   1.100 +	// (u, v) is the relative position of pt within the 4pix square
   1.101 +	// so if it's exactly at the lower-right it'll be (0, 0), in the middle it'll be (0.5, 0.5) etc.
   1.102 +	float u = (pt.x - x0) * fxsz;
   1.103 +	float v = (pt.y - y0) * fysz;
   1.104 +
   1.105 +	// we need integer pixel coordinates to use lookup()
   1.106 +	int px = (int)floorx;
   1.107 +	int py = (int)floory;
   1.108 +
   1.109 +	// the heights at the 4 corners
   1.110 +	float h00 = lookup(px, py);
   1.111 +	float h10 = lookup(px + 1, py);
   1.112 +	float h01 = lookup(px, py + 1);
   1.113 +	float h11 = lookup(px + 1, py + 1);
   1.114 +
   1.115 +	// first interpolate along the top and bottom edges
   1.116 +	float top = lerp(h01, h11, u);
   1.117 +	float bot = lerp(h00, h10, u);
   1.118 +	// then vertically between the top/bot values
   1.119 +	return lerp(bot, top, v);
   1.120 +}
   1.121 +
   1.122 +Vector3 Terrain::get_normal(const Vector2 &pt) const
   1.123 +{
   1.124 +	// compute partial derivatives (slopes) in respect to X and Y using central differences
   1.125 +	float dfdx = get_height(Vector2(pt.x + dx, pt.y)) - get_height(Vector2(pt.x - dx, pt.y));
   1.126 +	float dfdy = get_height(Vector2(pt.x, pt.y + dy)) - get_height(Vector2(pt.x, pt.y - dy));
   1.127 +
   1.128 +	/* we need to multiply with the final scale factors while constructing the tangent and
   1.129 +	 * bitangent vectors, otherwise the "aspect" of the normals will be wrong while changing
   1.130 +	 * the two scale factors independently.
   1.131 +	 */
   1.132 +	Vector3 tang = Vector3(2.0 * dx * scale, dfdx * height_scale, 0);
   1.133 +	Vector3 bitan = Vector3(0, dfdy * height_scale, 2.0 * dy * scale);
   1.134 +	return cross_product(bitan, tang).normalized();
   1.135 +}
   1.136 +
   1.137 +float Terrain::get_world_height(const Vector2 &pt) const
   1.138 +{
   1.139 +	return get_height(world_to_uv(pt)) * height_scale;
   1.140 +}
   1.141 +
   1.142 +/* step over pixel distances until we cross the terrain, then determine the
   1.143 + * point of intersection
   1.144 + */
   1.145 +bool Terrain::intersect(const Ray &ray, float *dist) const
   1.146 +{
   1.147 +	// axis-aligned bounding box of the terrain
   1.148 +	AABox aabb;
   1.149 +	aabb.min = Vector3(-0.5 * scale, 0, -0.5 * scale);
   1.150 +	aabb.max = Vector3(0.5 * scale, height_scale, 0.5 * scale);
   1.151 +
   1.152 +	Ray mray;
   1.153 +	mray.origin = ray.origin;
   1.154 +	// find a reasonable step size
   1.155 +	float pixsz = dx > dy ? dx : dy;
   1.156 +	float raysz = ray.dir.length();
   1.157 +	mray.dir = ray.dir / raysz * pixsz;
   1.158 +
   1.159 +	float aabb_dist = 0.0;
   1.160 +
   1.161 +	if(!aabb.contains(mray.origin)) {
   1.162 +		/* if we're not in the aabb, calculate the intersection point of the ray
   1.163 +		 * with the aabb, to start ray-marching from that point.
   1.164 +		 */
   1.165 +		HitPoint hit;
   1.166 +		if(!aabb.intersect(ray, &hit)) {
   1.167 +			return false;	// the ray misses the terrain completely
   1.168 +		}
   1.169 +		aabb_dist = hit.dist;
   1.170 +		mray.origin += ray.dir * aabb_dist;
   1.171 +	}
   1.172 +
   1.173 +
   1.174 +	// ray-march over the terrain until we find an intersection or leave the box
   1.175 +	int num_steps = 0;
   1.176 +	do {
   1.177 +		float mdist;
   1.178 +		if(intersect_micro(mray, &mdist, 0.01)) {
   1.179 +			/* to calculate the distance travelled along the original ray we need to
   1.180 +			 * first add the distance to the aabb where we started ray marching. Then
   1.181 +			 * add the fraction of the last step to the number of micro-steps up to
   1.182 +			 * this point, which have to be converted to the original ray scale by
   1.183 +			 * multiplying by the magnitude of the micro-ray direction, over the
   1.184 +			 * magnitude of the original ray direction.
   1.185 +			 */
   1.186 +			*dist = aabb_dist + (mdist + (float)num_steps) * pixsz / raysz;
   1.187 +			return true;
   1.188 +		}
   1.189 +
   1.190 +		mray.origin += mray.dir;
   1.191 +		num_steps++;
   1.192 +	} while(aabb.contains(mray.origin));
   1.193 +
   1.194 +	return false;
   1.195 +}
   1.196 +
   1.197 +bool Terrain::intersect_micro(const Ray &ray, float *dist, float thres) const
   1.198 +{
   1.199 +	Vector3 start = ray.origin;
   1.200 +	Vector3 end = ray.origin + ray.dir;
   1.201 +
   1.202 +	float hstart = get_world_height(Vector2(start.x, start.z));
   1.203 +	float hend = get_world_height(Vector2(end.x, end.z));
   1.204 +
   1.205 +	/* if both the start and the end of the ray are above or below the heightfield
   1.206 +	 * then there's no possible intersection in micro-scales.
   1.207 +	 */
   1.208 +	if(start.y > hstart && end.y > hend) {
   1.209 +		return false;	// all above
   1.210 +	}
   1.211 +	if(start.y < hstart && end.y < hend) {
   1.212 +		return false;	// all below
   1.213 +	}
   1.214 +
   1.215 +	/* otherwise (one on one side and the other on the other side), we're straddling
   1.216 +	 * the heightfield. So, find the mid-point and binary-search until we find an
   1.217 +	 * estimated intersection point within the limits defined by "thres".
   1.218 +	 */
   1.219 +	Vector3 mid = lerp(start, end, 0.5);
   1.220 +	float hmid = get_world_height(Vector2(mid.x, mid.z));
   1.221 +	float dh = mid.y - hmid;
   1.222 +
   1.223 +	int iter = 0;	// iter is there to avoid infinite loops in marginal cases
   1.224 +	while(fabs(dh) > thres && iter < 100) {
   1.225 +		if(dh > 0) {	// mid is above the surface
   1.226 +			start = mid;
   1.227 +		} else {		// mid is below the surface
   1.228 +			end = mid;
   1.229 +		}
   1.230 +		mid = lerp(start, end, 0.5);
   1.231 +		hmid = get_world_height(Vector2(mid.x, mid.z));
   1.232 +		dh = mid.y - hmid;
   1.233 +
   1.234 +		iter++;
   1.235 +	}
   1.236 +
   1.237 +	// found the intersection point, calculate the parametric distance and return true
   1.238 +	*dist = (mid - ray.origin).length() / ray.dir.length();
   1.239 +	return true;
   1.240 +}
   1.241 +
   1.242 +void Terrain::draw(int usub, int vsub) const
   1.243 +{
   1.244 +	if(usub < 1) usub = xsz - 1;
   1.245 +	if(vsub < 1) vsub = ysz - 1;
   1.246 +
   1.247 +	if(dlist && dlist_sub[0] == usub && dlist_sub[1] == vsub) {
   1.248 +		glCallList(dlist);
   1.249 +		return;
   1.250 +	}
   1.251 +
   1.252 +	if(!dlist) {
   1.253 +		dlist = glGenLists(2);
   1.254 +		dlist_norm = dlist + 1;
   1.255 +	}
   1.256 +	glNewList(dlist, GL_COMPILE_AND_EXECUTE);
   1.257 +	dlist_sub[0] = usub;
   1.258 +	dlist_sub[1] = vsub;
   1.259 +
   1.260 +	float du = 1.0 / (float)usub;
   1.261 +	float dv = 1.0 / (float)vsub;
   1.262 +
   1.263 +	glBegin(GL_QUADS);
   1.264 +
   1.265 +	float v = 0.0;
   1.266 +	for(int i=0; i<vsub; i++) {
   1.267 +		float u = 0.0;
   1.268 +		for(int j=0; j<usub; j++) {
   1.269 +
   1.270 +			Vector2 uv(u, v);
   1.271 +			Vector3 norm = get_normal(uv);
   1.272 +			Vector2 wpt = uv_to_world(Vector2(u, v));
   1.273 +			glNormal3f(norm.x, norm.y, norm.z);
   1.274 +			glVertex3f(wpt.x, get_height(Vector2(u, v)) * height_scale, wpt.y);
   1.275 +
   1.276 +			uv = Vector2(u, v + dv);
   1.277 +			norm = get_normal(uv);
   1.278 +			wpt = uv_to_world(uv);
   1.279 +			glNormal3f(norm.x, norm.y, norm.z);
   1.280 +			glVertex3f(wpt.x, get_height(uv) * height_scale, wpt.y);
   1.281 +
   1.282 +			uv = Vector2(u + du, v + dv);
   1.283 +			norm = get_normal(uv);
   1.284 +			wpt = uv_to_world(uv);
   1.285 +			glNormal3f(norm.x, norm.y, norm.z);
   1.286 +			glVertex3f(wpt.x, get_height(uv) * height_scale, wpt.y);
   1.287 +
   1.288 +			uv = Vector2(u + du, v);
   1.289 +			norm = get_normal(uv);
   1.290 +			wpt = uv_to_world(uv);
   1.291 +			glNormal3f(norm.x, norm.y, norm.z);
   1.292 +			glVertex3f(wpt.x, get_height(uv) * height_scale, wpt.y);
   1.293 +
   1.294 +			u += du;
   1.295 +		}
   1.296 +		v += dv;
   1.297 +	}
   1.298 +	glEnd();
   1.299 +
   1.300 +	glEndList();
   1.301 +}
   1.302 +
   1.303 +void Terrain::draw_normals(float sz, int usub, int vsub) const
   1.304 +{
   1.305 +	if(usub < 1) usub = xsz - 1;
   1.306 +	if(vsub < 1) vsub = ysz - 1;
   1.307 +
   1.308 +	if(dlist_norm && dlist_norm_sub[0] == usub && dlist_norm_sub[1] == vsub && fabs(dlist_norm_sz - sz) < 0.0001) {
   1.309 +		glCallList(dlist_norm);
   1.310 +		return;
   1.311 +	}
   1.312 +
   1.313 +	if(!dlist_norm) {
   1.314 +		dlist = glGenLists(2);
   1.315 +		dlist_norm = dlist + 1;
   1.316 +	}
   1.317 +	glNewList(dlist_norm, GL_COMPILE_AND_EXECUTE);
   1.318 +	dlist_norm_sub[0] = usub;
   1.319 +	dlist_norm_sub[1] = vsub;
   1.320 +	dlist_norm_sz = sz;
   1.321 +
   1.322 +	float du = 1.0 / (float)usub;
   1.323 +	float dv = 1.0 / (float)vsub;
   1.324 +
   1.325 +	glBegin(GL_LINES);
   1.326 +
   1.327 +	float v = 0.0;
   1.328 +	for(int i=0; i<vsub + 1; i++) {
   1.329 +		float u = 0.0;
   1.330 +		for(int j=0; j<usub + 1; j++) {
   1.331 +			Vector2 uv(u, v);
   1.332 +			Vector2 wpt = uv_to_world(uv);
   1.333 +			Vector3 norm = get_normal(uv) * sz * height_scale;
   1.334 +
   1.335 +			Vector3 p(wpt.x, get_height(uv) * height_scale, wpt.y);
   1.336 +			glVertex3f(p.x, p.y, p.z);
   1.337 +			glVertex3f(p.x + norm.x, p.y + norm.y, p.z + norm.z);
   1.338 +
   1.339 +			u += du;
   1.340 +		}
   1.341 +		v += dv;
   1.342 +	}
   1.343 +	glEnd();
   1.344 +
   1.345 +	glEndList();
   1.346 +}
   1.347 +
   1.348 +void Terrain::invalidate_mesh()
   1.349 +{
   1.350 +	// this will force recreation of the display lists on the next draw call
   1.351 +	dlist_sub[0] = dlist_sub[1] = dlist_norm_sub[0] = dlist_norm_sub[1] = -1;
   1.352 +}