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