medfall

A super great game engine
Log | Files | Refs

heightmap.cc (7779B)


      1 #include <stdio.h>
      2 
      3 #include "intrinsics.h"
      4 #include "linear_algebra.h"
      5 #include "aabb.h"
      6 #include "heightmap.h"
      7 #include "log.h"
      8 #include "memory_arena.h"
      9 #include "decompress_bc.h"
     10 
     11 static float bilerp( v3 p1, v3 p2, v3 p3, v3 p4, v2 p ) {
     12 	const float tx = ( p.x - p1.x ) / ( p2.x - p1.x );
     13 
     14 	const float mx1 = lerp( p1.z, tx, p2.z );
     15 	const float mx2 = lerp( p3.z, tx, p4.z );
     16 
     17 	const float ty = ( p.y - p1.y ) / ( p3.y - p1.y );
     18 
     19 	return lerp( mx1, ty, mx2 );
     20 }
     21 
     22 float heightmap_height( const array2d< u16 > hm, u32 x, u32 y ) {
     23 	return hm.try_get( x, y, 0 ) / 256.0f;
     24 }
     25 
     26 v3 heightmap_point( const array2d< u16 > hm, u32 x, u32 y ) {
     27 	return v3( checked_cast< float >( x ), checked_cast< float >( y ), heightmap_height( hm, x, y ) );
     28 }
     29 
     30 static bool ray_vs_aabb( const MinMax & aabb, const Ray3 & ray, float * t ) {
     31 	if( aabb.contains( ray.origin ) ) {
     32 		*t = 0.0f;
     33 		return true;
     34 	}
     35 
     36 	float tx1 = ( aabb.mins.x - ray.origin.x ) * ray.inv_dir.x;
     37 	float tx2 = ( aabb.maxs.x - ray.origin.x ) * ray.inv_dir.x;
     38 
     39 	float tmin = min( tx1, tx2 );
     40 	float tmax = max( tx1, tx2 );
     41 
     42 	float ty1 = ( aabb.mins.y - ray.origin.y ) * ray.inv_dir.y;
     43 	float ty2 = ( aabb.maxs.y - ray.origin.y ) * ray.inv_dir.y;
     44 
     45 	tmin = max( tmin, min( ty1, ty2 ) );
     46 	tmax = min( tmax, max( ty1, ty2 ) );
     47 
     48 	float tz1 = ( aabb.mins.z - ray.origin.z ) * ray.inv_dir.z;
     49 	float tz2 = ( aabb.maxs.z - ray.origin.z ) * ray.inv_dir.z;
     50 
     51 	tmin = max( tmin, min( tz1, tz2 ) );
     52 	tmax = min( tmax, max( tz1, tz2 ) );
     53 
     54 	if( tmax >= tmin ) {
     55 		*t = tmin;
     56 		return true;
     57 	}
     58 	return false;
     59 }
     60 
     61 bool ray_vs_sphere( v3 ray_origin, v3 ray_dir, v3 sphere_origin, float sphere_radius, float * t ) {
     62 	v3 m = ray_origin - sphere_origin;
     63 
     64 	float b = dot( m, ray_dir );
     65 	float c = dot( m, m ) - sphere_radius * sphere_radius;
     66 
     67 	float determinant = b * b - c;
     68 	if( determinant < 0 ) {
     69 		return false;
     70 	}
     71 
     72 	*t = -b - sqrtf( determinant );
     73 	return *t >= 0;
     74 }
     75 
     76 bool ray_vs_plane( v3 ray_origin, v3 ray_dir, v3 plane_normal, float plane_distance, float * t ) {
     77 	const float epsilon = 0.001f;
     78 	float denom = dot( ray_dir, plane_normal );
     79 	if( fabsf( denom ) < epsilon ) {
     80 		return false;
     81 	}
     82 
     83 	*t = ( plane_distance - dot( ray_origin, plane_normal ) ) / denom;
     84 	return *t >= 0;
     85 }
     86 
     87 void barycentric( v3 a, v3 b, v3 c, v3 p, float * u, float * v, float * w ) {
     88 	v3 ab = b - a;
     89 	v3 ac = c - a;
     90 	v3 ap = p - a;
     91 
     92 	float dabab = dot( ab, ab );
     93 	float dabac = dot( ab, ac );
     94 	float dacac = dot( ac, ac );
     95 	float dapab = dot( ap, ab );
     96 	float dapac = dot( ap, ac );
     97 	float denom = dabab * dacac - dabac * dabac;
     98 
     99 	*v = ( dacac * dapab - dabac * dapac ) / denom;
    100 	*w = ( dabab * dapac - dabac * dapab ) / denom;
    101 	*u = 1.0f - *v - *w;
    102 }
    103 
    104 bool point_in_triangle( v3 p, v3 a, v3 b, v3 c ) {
    105 	float u, v, w;
    106 	barycentric( a, b, c, p, &u, &v, &w );
    107 
    108 	const float epsilon = 0.0001f;
    109 	if( u < -epsilon || u > 1.0f + epsilon ) return false;
    110 	if( v < -epsilon || v > 1.0f + epsilon ) return false;
    111 	if( w < -epsilon || w > 1.0f + epsilon ) return false;
    112 
    113 	return true;
    114 }
    115 
    116 // TODO: not robust w/ thin triangles
    117 v3 triangle_normal( v3 a, v3 b, v3 c ) {
    118 	return normalize( cross( b - a, c - a ) );
    119 }
    120 
    121 bool ray_vs_triangle( v3 ray_origin, v3 ray_dir, v3 p0, v3 p1, v3 p2, float * t ) {
    122 	v3 e01 = p1 - p0;
    123 	v3 e02 = p2 - p0;
    124 	v3 p = cross( ray_dir, e02 );
    125 	float det = dot( e01, p );
    126 
    127 	const float epsilon = 0.001f;
    128 	if( fabsf( det ) < epsilon ) return false;
    129 
    130 	float inv_det = 1.0f / det;
    131 
    132 	v3 d = ray_origin - p0;
    133 	float u = dot( d, p ) * inv_det;
    134 	if( u < 0 || u > 1 ) return false;
    135 
    136 	v3 q = cross( d, e01 );
    137 	float v = dot( ray_dir, q ) * inv_det;
    138 	if( v < 0 || u + v > 1 ) return false;
    139 
    140 	*t = dot( e02, q ) * inv_det;
    141 	return *t >= 0;
    142 }
    143 
    144 static size_t child_idx( size_t node_idx, size_t quadrant ) {
    145 	return node_idx * 4 + quadrant + 1;
    146 }
    147 
    148 static void sort4( int * indices, float * elems ) {
    149 	for( int i = 0; i < 4; i++ )
    150 		indices[ i ] = i;
    151 
    152 	// sort a/b and c/d
    153 	if( elems[ indices[ 0 ] ] > elems[ indices[ 1 ] ] )
    154 		swap2( indices[ 0 ], indices[ 1 ] );
    155 	if( elems[ indices[ 2 ] ] > elems[ indices[ 3 ] ] )
    156 		swap2( indices[ 2 ], indices[ 3 ] );
    157 
    158 	// find lowest and highest
    159 	if( elems[ indices[ 0 ] ] > elems[ indices[ 2 ] ] )
    160 		swap2( indices[ 0 ], indices[ 2 ] );
    161 	if( elems[ indices[ 1 ] ] > elems[ indices[ 3 ] ] )
    162 		swap2( indices[ 1 ], indices[ 3 ] );
    163 
    164 	// sort middle two
    165 	if( elems[ indices[ 1 ] ] > elems[ indices[ 2 ] ] )
    166 		swap2( indices[ 1 ], indices[ 2 ] );
    167 }
    168 
    169 static bool ray_vs_quadtree_node(
    170 	const QuadTree & qt, size_t node_idx, MinMaxu32 aabb,
    171 	const Ray3 & ray, float * t, v3 * xnormal, float max_dist
    172 ) {
    173 	if( aabb.maxs.x - aabb.mins.x == 2 || aabb.maxs.y - aabb.mins.y == 2 ) {
    174 		bool hit = false;
    175 
    176 		for( size_t i = 0; i < 2; i++ ) {
    177 			for( size_t j = 0; j < 2; j++ ) {
    178 				u32 x = aabb.mins.x + j;
    179 				u32 y = aabb.mins.y + i;
    180 
    181 				v3 p0 = heightmap_point( qt.heightmap, x + 0, y + 0 );
    182 				v3 p1 = heightmap_point( qt.heightmap, x + 1, y + 0 );
    183 				v3 p2 = heightmap_point( qt.heightmap, x + 0, y + 1 );
    184 				v3 p3 = heightmap_point( qt.heightmap, x + 1, y + 1 );
    185 
    186 				// bottom right triangle
    187 				float t0;
    188 				v3 normal0 = triangle_normal( p0, p1, p3 );
    189 				bool hit0 = ray_vs_triangle( ray.origin, ray.dir, p0, p1, p3, &t0 );
    190 				if( hit0 ) {
    191 					if( !hit || t0 < *t ) {
    192 						*t = t0;
    193 						*xnormal = normal0;
    194 						hit = true;
    195 					}
    196 				}
    197 
    198 				// top left triangle
    199 				float t1;
    200 				v3 normal1 = triangle_normal( p0, p3, p2 );
    201 				bool hit1 = ray_vs_triangle( ray.origin, ray.dir, p0, p3, p2, &t1 );
    202 				if( hit1 ) {
    203 					if( !hit || t1 < *t ) {
    204 						*t = t1;
    205 						*xnormal = normal1;
    206 						hit = true;
    207 					}
    208 				}
    209 			}
    210 		}
    211 
    212 		return hit;
    213 	}
    214 
    215 	MinMaxu32 aabbs[ 4 ];
    216 	float ts[ 4 ];
    217 	for( int i = 0; i < 4; i++ ) {
    218 		QuadTreeNode child = qt.nodes[ child_idx( node_idx, i ) ];
    219 		aabbs[ i ] = aabb.quadrant( i ).clamp_z( child.min_z, child.max_z );
    220 		ts[ i ] = FLT_MAX;
    221 		ray_vs_aabb( MinMax( aabbs[ i ] ), ray, &ts[ i ] );
    222 	}
    223 
    224 	int sorted[ 4 ];
    225 	sort4( sorted, ts );
    226 
    227 	for( int i = 0; i < 4; i++ ) {
    228 		int idx = sorted[ i ];
    229 		if( ts[ idx ] >= max_dist )
    230 			break;
    231 
    232 		if( ray_vs_quadtree_node( qt, child_idx( node_idx, idx ), aabbs[ idx ], ray, t, xnormal, max_dist ) ) {
    233 			return true;
    234 		}
    235 	}
    236 
    237 	return false;
    238 }
    239 
    240 bool ray_vs_quadtree( const QuadTree & qt, const Ray3 & ray, float * t, v3 * xnormal, float max_dist ) {
    241 	v3u32 mins = v3u32( 0, 0, qt.nodes[ 0 ].min_z );
    242 	v3u32 maxs = v3u32( qt.heightmap.w, qt.heightmap.h, qt.nodes[ 0 ].max_z );
    243 	MinMaxu32 aabb( mins, maxs );
    244 
    245 	return ray_vs_quadtree_node( qt, 0, aabb, ray, t, xnormal, max_dist );
    246 }
    247 
    248 bool ray_vs_terrain( const QuadTree & qt, const Ray3 & ray, float * t, v3 * xnormal, float max_dist ) {
    249 	float dont_care_t;
    250 	v3 dont_care_normal;
    251 
    252 	if( t == NULL )
    253 		t = &dont_care_t;
    254 	if( xnormal == NULL )
    255 		xnormal = &dont_care_normal;
    256 
    257 	bool hit_quadtree = ray_vs_quadtree( qt, ray, t, xnormal, max_dist );
    258 	if( hit_quadtree )
    259 		return true;
    260 
    261 	*xnormal = v3( 0, 0, 1 );
    262 	return ray_vs_plane( ray.origin, ray.dir, *xnormal, 0, t );
    263 }
    264 
    265 bool segment_vs_terrain( const QuadTree & qt, const v3 & start, const v3 & end, float * t, v3 * xnormal ) {
    266 	if( t != NULL )
    267 		*t = 1.0f;
    268 
    269 	if( start == end )
    270 		return false;
    271 
    272 	Ray3 ray( start, normalize( end - start ) );
    273 	float l = length( end - start );
    274 
    275 	float absolute_t;
    276 	bool ok = ray_vs_terrain( qt, ray, &absolute_t, xnormal, l );
    277 	if( !ok || absolute_t >= l )
    278 		return false;
    279 
    280 	if( t != NULL )
    281 		*t = absolute_t / l;
    282 	return true;
    283 }
    284 
    285 void bc5_to_heightmap( MemoryArena * arena, array2d< u16 > hm, u8 * bc5 ) {
    286 	MEMARENA_SCOPED_CHECKPOINT( arena );
    287 	array2d< v2 > rg = alloc_array2d< v2 >( arena, hm.w, hm.h );
    288 	decompress_bc5( rg, bc5 );
    289 
    290 	for( u32 y = 0; y < rg.h; y++ ) {
    291 		for( u32 x = 0; x < rg.w; x++ ) {
    292 			float r = rg( x, y ).x;
    293 			float g = rg( x, y ).y;
    294 			hm( x, y ) = u16( ( r * 256.0f + g ) * 255.0f );
    295 		}
    296 	}
    297 }