nuclear@1: #include nuclear@1: #include "terrain.h" nuclear@1: #include "imago2.h" nuclear@1: nuclear@1: using namespace goatgfx; nuclear@1: nuclear@1: Terrain::Terrain() nuclear@1: { nuclear@1: xsz = ysz = 0; nuclear@1: dx = dy = 0; nuclear@1: hmap = 0; nuclear@1: scale = height_scale = 1; nuclear@1: nuclear@1: dlist = dlist_norm = 0; nuclear@1: nuclear@1: dlist_sub[0] = dlist_sub[1] = dlist_norm_sub[0] = dlist_norm_sub[1] = -1; nuclear@1: dlist_norm_sz = -1.0; nuclear@1: } nuclear@1: nuclear@1: Terrain::~Terrain() nuclear@1: { nuclear@1: if(hmap) { nuclear@1: img_free_pixels(hmap); nuclear@1: } nuclear@1: nuclear@1: glDeleteLists(dlist, 2); nuclear@1: } nuclear@1: nuclear@1: void Terrain::set_scale(float s) nuclear@1: { nuclear@1: if(s != scale) { nuclear@1: scale = s; nuclear@1: invalidate_mesh(); nuclear@1: } nuclear@1: } nuclear@1: nuclear@1: float Terrain::get_scale() const nuclear@1: { nuclear@1: return scale; nuclear@1: } nuclear@1: nuclear@1: void Terrain::set_height_scale(float s) nuclear@1: { nuclear@1: if(s != height_scale) { nuclear@1: height_scale = s; nuclear@1: invalidate_mesh(); nuclear@1: } nuclear@1: } nuclear@1: nuclear@1: float Terrain::get_height_scale() const nuclear@1: { nuclear@1: return height_scale; nuclear@1: } nuclear@1: nuclear@1: bool Terrain::load(const char *fname) nuclear@1: { nuclear@1: if(!(hmap = (float*)img_load_pixels(fname, &xsz, &ysz, IMG_FMT_GREYF))) { nuclear@1: fprintf(stderr, "failed to load image: %s\n", fname); nuclear@1: return false; nuclear@1: } nuclear@1: nuclear@1: dx = 1.0 / (float)(xsz - 1); nuclear@1: dy = 1.0 / (float)(ysz - 1); nuclear@1: return true; nuclear@1: } nuclear@1: nuclear@1: Vector2 Terrain::world_to_uv(const Vector2 &pt) const nuclear@1: { nuclear@1: return pt / scale + 0.5; nuclear@1: } nuclear@1: nuclear@1: Vector2 Terrain::uv_to_world(const Vector2 &pt) const nuclear@1: { nuclear@1: return (pt - 0.5) * scale; nuclear@1: } nuclear@1: nuclear@1: float Terrain::lookup(int x, int y) const nuclear@1: { nuclear@1: int px = x < 0 ? 0 : (x >= xsz ? x = xsz - 1 : x); nuclear@1: int py = y < 0 ? 0 : (y >= ysz ? y = ysz - 1 : y); nuclear@1: return hmap[py * xsz + px]; nuclear@1: } nuclear@1: nuclear@1: float Terrain::get_height(const Vector2 &pt) const nuclear@1: { nuclear@1: float fxsz = (float)xsz; nuclear@1: float fysz = (float)ysz; nuclear@1: nuclear@1: // the floor of x * xsz is the pixel x coordinate (e.g. [0.0, 1.0, 2.0, ... 127.0] for 128x128) nuclear@1: float floorx = floor(pt.x * fxsz); nuclear@1: float floory = floor(pt.y * fysz); nuclear@1: nuclear@1: // (x0, y0) is the lower left pixel in normalized coords of the four pixels surrounding pt nuclear@1: float x0 = floorx / fxsz; nuclear@1: float y0 = floory / fysz; nuclear@1: nuclear@1: // (u, v) is the relative position of pt within the 4pix square nuclear@1: // 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. nuclear@1: float u = (pt.x - x0) * fxsz; nuclear@1: float v = (pt.y - y0) * fysz; nuclear@1: nuclear@1: // we need integer pixel coordinates to use lookup() nuclear@1: int px = (int)floorx; nuclear@1: int py = (int)floory; nuclear@1: nuclear@1: // the heights at the 4 corners nuclear@1: float h00 = lookup(px, py); nuclear@1: float h10 = lookup(px + 1, py); nuclear@1: float h01 = lookup(px, py + 1); nuclear@1: float h11 = lookup(px + 1, py + 1); nuclear@1: nuclear@1: // first interpolate along the top and bottom edges nuclear@1: float top = lerp(h01, h11, u); nuclear@1: float bot = lerp(h00, h10, u); nuclear@1: // then vertically between the top/bot values nuclear@1: return lerp(bot, top, v); nuclear@1: } nuclear@1: nuclear@1: Vector3 Terrain::get_normal(const Vector2 &pt) const nuclear@1: { nuclear@1: // compute partial derivatives (slopes) in respect to X and Y using central differences nuclear@1: float dfdx = get_height(Vector2(pt.x + dx, pt.y)) - get_height(Vector2(pt.x - dx, pt.y)); nuclear@1: float dfdy = get_height(Vector2(pt.x, pt.y + dy)) - get_height(Vector2(pt.x, pt.y - dy)); nuclear@1: nuclear@1: /* we need to multiply with the final scale factors while constructing the tangent and nuclear@1: * bitangent vectors, otherwise the "aspect" of the normals will be wrong while changing nuclear@1: * the two scale factors independently. nuclear@1: */ nuclear@1: Vector3 tang = Vector3(2.0 * dx * scale, dfdx * height_scale, 0); nuclear@1: Vector3 bitan = Vector3(0, dfdy * height_scale, 2.0 * dy * scale); nuclear@1: return cross_product(bitan, tang).normalized(); nuclear@1: } nuclear@1: nuclear@1: float Terrain::get_world_height(const Vector2 &pt) const nuclear@1: { nuclear@1: return get_height(world_to_uv(pt)) * height_scale; nuclear@1: } nuclear@1: nuclear@1: /* step over pixel distances until we cross the terrain, then determine the nuclear@1: * point of intersection nuclear@1: */ nuclear@1: bool Terrain::intersect(const Ray &ray, float *dist) const nuclear@1: { nuclear@1: // axis-aligned bounding box of the terrain nuclear@1: AABox aabb; nuclear@1: aabb.min = Vector3(-0.5 * scale, 0, -0.5 * scale); nuclear@1: aabb.max = Vector3(0.5 * scale, height_scale, 0.5 * scale); nuclear@1: nuclear@1: Ray mray; nuclear@1: mray.origin = ray.origin; nuclear@1: // find a reasonable step size nuclear@1: float pixsz = dx > dy ? dx : dy; nuclear@1: float raysz = ray.dir.length(); nuclear@1: mray.dir = ray.dir / raysz * pixsz; nuclear@1: nuclear@1: float aabb_dist = 0.0; nuclear@1: nuclear@1: if(!aabb.contains(mray.origin)) { nuclear@1: /* if we're not in the aabb, calculate the intersection point of the ray nuclear@1: * with the aabb, to start ray-marching from that point. nuclear@1: */ nuclear@1: HitPoint hit; nuclear@1: if(!aabb.intersect(ray, &hit)) { nuclear@1: return false; // the ray misses the terrain completely nuclear@1: } nuclear@1: aabb_dist = hit.dist; nuclear@1: mray.origin += ray.dir * aabb_dist; nuclear@1: } nuclear@1: nuclear@1: nuclear@1: // ray-march over the terrain until we find an intersection or leave the box nuclear@1: int num_steps = 0; nuclear@1: do { nuclear@1: float mdist; nuclear@1: if(intersect_micro(mray, &mdist, 0.01)) { nuclear@1: /* to calculate the distance travelled along the original ray we need to nuclear@1: * first add the distance to the aabb where we started ray marching. Then nuclear@1: * add the fraction of the last step to the number of micro-steps up to nuclear@1: * this point, which have to be converted to the original ray scale by nuclear@1: * multiplying by the magnitude of the micro-ray direction, over the nuclear@1: * magnitude of the original ray direction. nuclear@1: */ nuclear@1: *dist = aabb_dist + (mdist + (float)num_steps) * pixsz / raysz; nuclear@1: return true; nuclear@1: } nuclear@1: nuclear@1: mray.origin += mray.dir; nuclear@1: num_steps++; nuclear@1: } while(aabb.contains(mray.origin)); nuclear@1: nuclear@1: return false; nuclear@1: } nuclear@1: nuclear@1: bool Terrain::intersect_micro(const Ray &ray, float *dist, float thres) const nuclear@1: { nuclear@1: Vector3 start = ray.origin; nuclear@1: Vector3 end = ray.origin + ray.dir; nuclear@1: nuclear@1: float hstart = get_world_height(Vector2(start.x, start.z)); nuclear@1: float hend = get_world_height(Vector2(end.x, end.z)); nuclear@1: nuclear@1: /* if both the start and the end of the ray are above or below the heightfield nuclear@1: * then there's no possible intersection in micro-scales. nuclear@1: */ nuclear@1: if(start.y > hstart && end.y > hend) { nuclear@1: return false; // all above nuclear@1: } nuclear@1: if(start.y < hstart && end.y < hend) { nuclear@1: return false; // all below nuclear@1: } nuclear@1: nuclear@1: /* otherwise (one on one side and the other on the other side), we're straddling nuclear@1: * the heightfield. So, find the mid-point and binary-search until we find an nuclear@1: * estimated intersection point within the limits defined by "thres". nuclear@1: */ nuclear@1: Vector3 mid = lerp(start, end, 0.5); nuclear@1: float hmid = get_world_height(Vector2(mid.x, mid.z)); nuclear@1: float dh = mid.y - hmid; nuclear@1: nuclear@1: int iter = 0; // iter is there to avoid infinite loops in marginal cases nuclear@1: while(fabs(dh) > thres && iter < 100) { nuclear@1: if(dh > 0) { // mid is above the surface nuclear@1: start = mid; nuclear@1: } else { // mid is below the surface nuclear@1: end = mid; nuclear@1: } nuclear@1: mid = lerp(start, end, 0.5); nuclear@1: hmid = get_world_height(Vector2(mid.x, mid.z)); nuclear@1: dh = mid.y - hmid; nuclear@1: nuclear@1: iter++; nuclear@1: } nuclear@1: nuclear@1: // found the intersection point, calculate the parametric distance and return true nuclear@1: *dist = (mid - ray.origin).length() / ray.dir.length(); nuclear@1: return true; nuclear@1: } nuclear@1: nuclear@1: void Terrain::draw(int usub, int vsub) const nuclear@1: { nuclear@1: if(usub < 1) usub = xsz - 1; nuclear@1: if(vsub < 1) vsub = ysz - 1; nuclear@1: nuclear@1: if(dlist && dlist_sub[0] == usub && dlist_sub[1] == vsub) { nuclear@1: glCallList(dlist); nuclear@1: return; nuclear@1: } nuclear@1: nuclear@1: if(!dlist) { nuclear@1: dlist = glGenLists(2); nuclear@1: dlist_norm = dlist + 1; nuclear@1: } nuclear@1: glNewList(dlist, GL_COMPILE_AND_EXECUTE); nuclear@1: dlist_sub[0] = usub; nuclear@1: dlist_sub[1] = vsub; nuclear@1: nuclear@1: float du = 1.0 / (float)usub; nuclear@1: float dv = 1.0 / (float)vsub; nuclear@1: nuclear@1: glBegin(GL_QUADS); nuclear@1: nuclear@1: float v = 0.0; nuclear@1: for(int i=0; i