symmetry

annotate src/terrain.cc @ 1:46fe847bba08

find changesets by author, revision, files, or words in the commit message
using goat3dgfx
author John Tsiombikas <nuclear@member.fsf.org>
date Tue, 25 Feb 2014 23:47:20 +0200
parents
children
rev   line source
nuclear@1 1 #include <goat3dgfx/goat3dgfx.h>
nuclear@1 2 #include "terrain.h"
nuclear@1 3 #include "imago2.h"
nuclear@1 4
nuclear@1 5 using namespace goatgfx;
nuclear@1 6
nuclear@1 7 Terrain::Terrain()
nuclear@1 8 {
nuclear@1 9 xsz = ysz = 0;
nuclear@1 10 dx = dy = 0;
nuclear@1 11 hmap = 0;
nuclear@1 12 scale = height_scale = 1;
nuclear@1 13
nuclear@1 14 dlist = dlist_norm = 0;
nuclear@1 15
nuclear@1 16 dlist_sub[0] = dlist_sub[1] = dlist_norm_sub[0] = dlist_norm_sub[1] = -1;
nuclear@1 17 dlist_norm_sz = -1.0;
nuclear@1 18 }
nuclear@1 19
nuclear@1 20 Terrain::~Terrain()
nuclear@1 21 {
nuclear@1 22 if(hmap) {
nuclear@1 23 img_free_pixels(hmap);
nuclear@1 24 }
nuclear@1 25
nuclear@1 26 glDeleteLists(dlist, 2);
nuclear@1 27 }
nuclear@1 28
nuclear@1 29 void Terrain::set_scale(float s)
nuclear@1 30 {
nuclear@1 31 if(s != scale) {
nuclear@1 32 scale = s;
nuclear@1 33 invalidate_mesh();
nuclear@1 34 }
nuclear@1 35 }
nuclear@1 36
nuclear@1 37 float Terrain::get_scale() const
nuclear@1 38 {
nuclear@1 39 return scale;
nuclear@1 40 }
nuclear@1 41
nuclear@1 42 void Terrain::set_height_scale(float s)
nuclear@1 43 {
nuclear@1 44 if(s != height_scale) {
nuclear@1 45 height_scale = s;
nuclear@1 46 invalidate_mesh();
nuclear@1 47 }
nuclear@1 48 }
nuclear@1 49
nuclear@1 50 float Terrain::get_height_scale() const
nuclear@1 51 {
nuclear@1 52 return height_scale;
nuclear@1 53 }
nuclear@1 54
nuclear@1 55 bool Terrain::load(const char *fname)
nuclear@1 56 {
nuclear@1 57 if(!(hmap = (float*)img_load_pixels(fname, &xsz, &ysz, IMG_FMT_GREYF))) {
nuclear@1 58 fprintf(stderr, "failed to load image: %s\n", fname);
nuclear@1 59 return false;
nuclear@1 60 }
nuclear@1 61
nuclear@1 62 dx = 1.0 / (float)(xsz - 1);
nuclear@1 63 dy = 1.0 / (float)(ysz - 1);
nuclear@1 64 return true;
nuclear@1 65 }
nuclear@1 66
nuclear@1 67 Vector2 Terrain::world_to_uv(const Vector2 &pt) const
nuclear@1 68 {
nuclear@1 69 return pt / scale + 0.5;
nuclear@1 70 }
nuclear@1 71
nuclear@1 72 Vector2 Terrain::uv_to_world(const Vector2 &pt) const
nuclear@1 73 {
nuclear@1 74 return (pt - 0.5) * scale;
nuclear@1 75 }
nuclear@1 76
nuclear@1 77 float Terrain::lookup(int x, int y) const
nuclear@1 78 {
nuclear@1 79 int px = x < 0 ? 0 : (x >= xsz ? x = xsz - 1 : x);
nuclear@1 80 int py = y < 0 ? 0 : (y >= ysz ? y = ysz - 1 : y);
nuclear@1 81 return hmap[py * xsz + px];
nuclear@1 82 }
nuclear@1 83
nuclear@1 84 float Terrain::get_height(const Vector2 &pt) const
nuclear@1 85 {
nuclear@1 86 float fxsz = (float)xsz;
nuclear@1 87 float fysz = (float)ysz;
nuclear@1 88
nuclear@1 89 // the floor of x * xsz is the pixel x coordinate (e.g. [0.0, 1.0, 2.0, ... 127.0] for 128x128)
nuclear@1 90 float floorx = floor(pt.x * fxsz);
nuclear@1 91 float floory = floor(pt.y * fysz);
nuclear@1 92
nuclear@1 93 // (x0, y0) is the lower left pixel in normalized coords of the four pixels surrounding pt
nuclear@1 94 float x0 = floorx / fxsz;
nuclear@1 95 float y0 = floory / fysz;
nuclear@1 96
nuclear@1 97 // (u, v) is the relative position of pt within the 4pix square
nuclear@1 98 // 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 99 float u = (pt.x - x0) * fxsz;
nuclear@1 100 float v = (pt.y - y0) * fysz;
nuclear@1 101
nuclear@1 102 // we need integer pixel coordinates to use lookup()
nuclear@1 103 int px = (int)floorx;
nuclear@1 104 int py = (int)floory;
nuclear@1 105
nuclear@1 106 // the heights at the 4 corners
nuclear@1 107 float h00 = lookup(px, py);
nuclear@1 108 float h10 = lookup(px + 1, py);
nuclear@1 109 float h01 = lookup(px, py + 1);
nuclear@1 110 float h11 = lookup(px + 1, py + 1);
nuclear@1 111
nuclear@1 112 // first interpolate along the top and bottom edges
nuclear@1 113 float top = lerp(h01, h11, u);
nuclear@1 114 float bot = lerp(h00, h10, u);
nuclear@1 115 // then vertically between the top/bot values
nuclear@1 116 return lerp(bot, top, v);
nuclear@1 117 }
nuclear@1 118
nuclear@1 119 Vector3 Terrain::get_normal(const Vector2 &pt) const
nuclear@1 120 {
nuclear@1 121 // compute partial derivatives (slopes) in respect to X and Y using central differences
nuclear@1 122 float dfdx = get_height(Vector2(pt.x + dx, pt.y)) - get_height(Vector2(pt.x - dx, pt.y));
nuclear@1 123 float dfdy = get_height(Vector2(pt.x, pt.y + dy)) - get_height(Vector2(pt.x, pt.y - dy));
nuclear@1 124
nuclear@1 125 /* we need to multiply with the final scale factors while constructing the tangent and
nuclear@1 126 * bitangent vectors, otherwise the "aspect" of the normals will be wrong while changing
nuclear@1 127 * the two scale factors independently.
nuclear@1 128 */
nuclear@1 129 Vector3 tang = Vector3(2.0 * dx * scale, dfdx * height_scale, 0);
nuclear@1 130 Vector3 bitan = Vector3(0, dfdy * height_scale, 2.0 * dy * scale);
nuclear@1 131 return cross_product(bitan, tang).normalized();
nuclear@1 132 }
nuclear@1 133
nuclear@1 134 float Terrain::get_world_height(const Vector2 &pt) const
nuclear@1 135 {
nuclear@1 136 return get_height(world_to_uv(pt)) * height_scale;
nuclear@1 137 }
nuclear@1 138
nuclear@1 139 /* step over pixel distances until we cross the terrain, then determine the
nuclear@1 140 * point of intersection
nuclear@1 141 */
nuclear@1 142 bool Terrain::intersect(const Ray &ray, float *dist) const
nuclear@1 143 {
nuclear@1 144 // axis-aligned bounding box of the terrain
nuclear@1 145 AABox aabb;
nuclear@1 146 aabb.min = Vector3(-0.5 * scale, 0, -0.5 * scale);
nuclear@1 147 aabb.max = Vector3(0.5 * scale, height_scale, 0.5 * scale);
nuclear@1 148
nuclear@1 149 Ray mray;
nuclear@1 150 mray.origin = ray.origin;
nuclear@1 151 // find a reasonable step size
nuclear@1 152 float pixsz = dx > dy ? dx : dy;
nuclear@1 153 float raysz = ray.dir.length();
nuclear@1 154 mray.dir = ray.dir / raysz * pixsz;
nuclear@1 155
nuclear@1 156 float aabb_dist = 0.0;
nuclear@1 157
nuclear@1 158 if(!aabb.contains(mray.origin)) {
nuclear@1 159 /* if we're not in the aabb, calculate the intersection point of the ray
nuclear@1 160 * with the aabb, to start ray-marching from that point.
nuclear@1 161 */
nuclear@1 162 HitPoint hit;
nuclear@1 163 if(!aabb.intersect(ray, &hit)) {
nuclear@1 164 return false; // the ray misses the terrain completely
nuclear@1 165 }
nuclear@1 166 aabb_dist = hit.dist;
nuclear@1 167 mray.origin += ray.dir * aabb_dist;
nuclear@1 168 }
nuclear@1 169
nuclear@1 170
nuclear@1 171 // ray-march over the terrain until we find an intersection or leave the box
nuclear@1 172 int num_steps = 0;
nuclear@1 173 do {
nuclear@1 174 float mdist;
nuclear@1 175 if(intersect_micro(mray, &mdist, 0.01)) {
nuclear@1 176 /* to calculate the distance travelled along the original ray we need to
nuclear@1 177 * first add the distance to the aabb where we started ray marching. Then
nuclear@1 178 * add the fraction of the last step to the number of micro-steps up to
nuclear@1 179 * this point, which have to be converted to the original ray scale by
nuclear@1 180 * multiplying by the magnitude of the micro-ray direction, over the
nuclear@1 181 * magnitude of the original ray direction.
nuclear@1 182 */
nuclear@1 183 *dist = aabb_dist + (mdist + (float)num_steps) * pixsz / raysz;
nuclear@1 184 return true;
nuclear@1 185 }
nuclear@1 186
nuclear@1 187 mray.origin += mray.dir;
nuclear@1 188 num_steps++;
nuclear@1 189 } while(aabb.contains(mray.origin));
nuclear@1 190
nuclear@1 191 return false;
nuclear@1 192 }
nuclear@1 193
nuclear@1 194 bool Terrain::intersect_micro(const Ray &ray, float *dist, float thres) const
nuclear@1 195 {
nuclear@1 196 Vector3 start = ray.origin;
nuclear@1 197 Vector3 end = ray.origin + ray.dir;
nuclear@1 198
nuclear@1 199 float hstart = get_world_height(Vector2(start.x, start.z));
nuclear@1 200 float hend = get_world_height(Vector2(end.x, end.z));
nuclear@1 201
nuclear@1 202 /* if both the start and the end of the ray are above or below the heightfield
nuclear@1 203 * then there's no possible intersection in micro-scales.
nuclear@1 204 */
nuclear@1 205 if(start.y > hstart && end.y > hend) {
nuclear@1 206 return false; // all above
nuclear@1 207 }
nuclear@1 208 if(start.y < hstart && end.y < hend) {
nuclear@1 209 return false; // all below
nuclear@1 210 }
nuclear@1 211
nuclear@1 212 /* otherwise (one on one side and the other on the other side), we're straddling
nuclear@1 213 * the heightfield. So, find the mid-point and binary-search until we find an
nuclear@1 214 * estimated intersection point within the limits defined by "thres".
nuclear@1 215 */
nuclear@1 216 Vector3 mid = lerp(start, end, 0.5);
nuclear@1 217 float hmid = get_world_height(Vector2(mid.x, mid.z));
nuclear@1 218 float dh = mid.y - hmid;
nuclear@1 219
nuclear@1 220 int iter = 0; // iter is there to avoid infinite loops in marginal cases
nuclear@1 221 while(fabs(dh) > thres && iter < 100) {
nuclear@1 222 if(dh > 0) { // mid is above the surface
nuclear@1 223 start = mid;
nuclear@1 224 } else { // mid is below the surface
nuclear@1 225 end = mid;
nuclear@1 226 }
nuclear@1 227 mid = lerp(start, end, 0.5);
nuclear@1 228 hmid = get_world_height(Vector2(mid.x, mid.z));
nuclear@1 229 dh = mid.y - hmid;
nuclear@1 230
nuclear@1 231 iter++;
nuclear@1 232 }
nuclear@1 233
nuclear@1 234 // found the intersection point, calculate the parametric distance and return true
nuclear@1 235 *dist = (mid - ray.origin).length() / ray.dir.length();
nuclear@1 236 return true;
nuclear@1 237 }
nuclear@1 238
nuclear@1 239 void Terrain::draw(int usub, int vsub) const
nuclear@1 240 {
nuclear@1 241 if(usub < 1) usub = xsz - 1;
nuclear@1 242 if(vsub < 1) vsub = ysz - 1;
nuclear@1 243
nuclear@1 244 if(dlist && dlist_sub[0] == usub && dlist_sub[1] == vsub) {
nuclear@1 245 glCallList(dlist);
nuclear@1 246 return;
nuclear@1 247 }
nuclear@1 248
nuclear@1 249 if(!dlist) {
nuclear@1 250 dlist = glGenLists(2);
nuclear@1 251 dlist_norm = dlist + 1;
nuclear@1 252 }
nuclear@1 253 glNewList(dlist, GL_COMPILE_AND_EXECUTE);
nuclear@1 254 dlist_sub[0] = usub;
nuclear@1 255 dlist_sub[1] = vsub;
nuclear@1 256
nuclear@1 257 float du = 1.0 / (float)usub;
nuclear@1 258 float dv = 1.0 / (float)vsub;
nuclear@1 259
nuclear@1 260 glBegin(GL_QUADS);
nuclear@1 261
nuclear@1 262 float v = 0.0;
nuclear@1 263 for(int i=0; i<vsub; i++) {
nuclear@1 264 float u = 0.0;
nuclear@1 265 for(int j=0; j<usub; j++) {
nuclear@1 266
nuclear@1 267 Vector2 uv(u, v);
nuclear@1 268 Vector3 norm = get_normal(uv);
nuclear@1 269 Vector2 wpt = uv_to_world(Vector2(u, v));
nuclear@1 270 glNormal3f(norm.x, norm.y, norm.z);
nuclear@1 271 glVertex3f(wpt.x, get_height(Vector2(u, v)) * height_scale, wpt.y);
nuclear@1 272
nuclear@1 273 uv = Vector2(u, v + dv);
nuclear@1 274 norm = get_normal(uv);
nuclear@1 275 wpt = uv_to_world(uv);
nuclear@1 276 glNormal3f(norm.x, norm.y, norm.z);
nuclear@1 277 glVertex3f(wpt.x, get_height(uv) * height_scale, wpt.y);
nuclear@1 278
nuclear@1 279 uv = Vector2(u + du, v + dv);
nuclear@1 280 norm = get_normal(uv);
nuclear@1 281 wpt = uv_to_world(uv);
nuclear@1 282 glNormal3f(norm.x, norm.y, norm.z);
nuclear@1 283 glVertex3f(wpt.x, get_height(uv) * height_scale, wpt.y);
nuclear@1 284
nuclear@1 285 uv = Vector2(u + du, v);
nuclear@1 286 norm = get_normal(uv);
nuclear@1 287 wpt = uv_to_world(uv);
nuclear@1 288 glNormal3f(norm.x, norm.y, norm.z);
nuclear@1 289 glVertex3f(wpt.x, get_height(uv) * height_scale, wpt.y);
nuclear@1 290
nuclear@1 291 u += du;
nuclear@1 292 }
nuclear@1 293 v += dv;
nuclear@1 294 }
nuclear@1 295 glEnd();
nuclear@1 296
nuclear@1 297 glEndList();
nuclear@1 298 }
nuclear@1 299
nuclear@1 300 void Terrain::draw_normals(float sz, int usub, int vsub) const
nuclear@1 301 {
nuclear@1 302 if(usub < 1) usub = xsz - 1;
nuclear@1 303 if(vsub < 1) vsub = ysz - 1;
nuclear@1 304
nuclear@1 305 if(dlist_norm && dlist_norm_sub[0] == usub && dlist_norm_sub[1] == vsub && fabs(dlist_norm_sz - sz) < 0.0001) {
nuclear@1 306 glCallList(dlist_norm);
nuclear@1 307 return;
nuclear@1 308 }
nuclear@1 309
nuclear@1 310 if(!dlist_norm) {
nuclear@1 311 dlist = glGenLists(2);
nuclear@1 312 dlist_norm = dlist + 1;
nuclear@1 313 }
nuclear@1 314 glNewList(dlist_norm, GL_COMPILE_AND_EXECUTE);
nuclear@1 315 dlist_norm_sub[0] = usub;
nuclear@1 316 dlist_norm_sub[1] = vsub;
nuclear@1 317 dlist_norm_sz = sz;
nuclear@1 318
nuclear@1 319 float du = 1.0 / (float)usub;
nuclear@1 320 float dv = 1.0 / (float)vsub;
nuclear@1 321
nuclear@1 322 glBegin(GL_LINES);
nuclear@1 323
nuclear@1 324 float v = 0.0;
nuclear@1 325 for(int i=0; i<vsub + 1; i++) {
nuclear@1 326 float u = 0.0;
nuclear@1 327 for(int j=0; j<usub + 1; j++) {
nuclear@1 328 Vector2 uv(u, v);
nuclear@1 329 Vector2 wpt = uv_to_world(uv);
nuclear@1 330 Vector3 norm = get_normal(uv) * sz * height_scale;
nuclear@1 331
nuclear@1 332 Vector3 p(wpt.x, get_height(uv) * height_scale, wpt.y);
nuclear@1 333 glVertex3f(p.x, p.y, p.z);
nuclear@1 334 glVertex3f(p.x + norm.x, p.y + norm.y, p.z + norm.z);
nuclear@1 335
nuclear@1 336 u += du;
nuclear@1 337 }
nuclear@1 338 v += dv;
nuclear@1 339 }
nuclear@1 340 glEnd();
nuclear@1 341
nuclear@1 342 glEndList();
nuclear@1 343 }
nuclear@1 344
nuclear@1 345 void Terrain::invalidate_mesh()
nuclear@1 346 {
nuclear@1 347 // this will force recreation of the display lists on the next draw call
nuclear@1 348 dlist_sub[0] = dlist_sub[1] = dlist_norm_sub[0] = dlist_norm_sub[1] = -1;
nuclear@1 349 }