medfall

heightmap.cc (7779B)

---

 #include <stdio.h>
 
 #include "intrinsics.h"
 #include "linear_algebra.h"
 #include "aabb.h"
 #include "heightmap.h"
 #include "log.h"
 #include "memory_arena.h"
 #include "decompress_bc.h"
 
 static float bilerp( v3 p1, v3 p2, v3 p3, v3 p4, v2 p ) {
 	const float tx = ( p.x - p1.x ) / ( p2.x - p1.x );
 
 	const float mx1 = lerp( p1.z, tx, p2.z );
 	const float mx2 = lerp( p3.z, tx, p4.z );
 
 	const float ty = ( p.y - p1.y ) / ( p3.y - p1.y );
 
 	return lerp( mx1, ty, mx2 );
 }
 
 float heightmap_height( const array2d< u16 > hm, u32 x, u32 y ) {
 	return hm.try_get( x, y, 0 ) / 256.0f;
 }
 
 v3 heightmap_point( const array2d< u16 > hm, u32 x, u32 y ) {
 	return v3( checked_cast< float >( x ), checked_cast< float >( y ), heightmap_height( hm, x, y ) );
 }
 
 static bool ray_vs_aabb( const MinMax & aabb, const Ray3 & ray, float * t ) {
 	if( aabb.contains( ray.origin ) ) {
 		*t = 0.0f;
 		return true;
 	}
 
 	float tx1 = ( aabb.mins.x - ray.origin.x ) * ray.inv_dir.x;
 	float tx2 = ( aabb.maxs.x - ray.origin.x ) * ray.inv_dir.x;
 
 	float tmin = min( tx1, tx2 );
 	float tmax = max( tx1, tx2 );
 
 	float ty1 = ( aabb.mins.y - ray.origin.y ) * ray.inv_dir.y;
 	float ty2 = ( aabb.maxs.y - ray.origin.y ) * ray.inv_dir.y;
 
 	tmin = max( tmin, min( ty1, ty2 ) );
 	tmax = min( tmax, max( ty1, ty2 ) );
 
 	float tz1 = ( aabb.mins.z - ray.origin.z ) * ray.inv_dir.z;
 	float tz2 = ( aabb.maxs.z - ray.origin.z ) * ray.inv_dir.z;
 
 	tmin = max( tmin, min( tz1, tz2 ) );
 	tmax = min( tmax, max( tz1, tz2 ) );
 
 	if( tmax >= tmin ) {
 		*t = tmin;
 		return true;
 	}
 	return false;
 }
 
 bool ray_vs_sphere( v3 ray_origin, v3 ray_dir, v3 sphere_origin, float sphere_radius, float * t ) {
 	v3 m = ray_origin - sphere_origin;
 
 	float b = dot( m, ray_dir );
 	float c = dot( m, m ) - sphere_radius * sphere_radius;
 
 	float determinant = b * b - c;
 	if( determinant < 0 ) {
 		return false;
 	}
 
 	*t = -b - sqrtf( determinant );
 	return *t >= 0;
 }
 
 bool ray_vs_plane( v3 ray_origin, v3 ray_dir, v3 plane_normal, float plane_distance, float * t ) {
 	const float epsilon = 0.001f;
 	float denom = dot( ray_dir, plane_normal );
 	if( fabsf( denom ) < epsilon ) {
 		return false;
 	}
 
 	*t = ( plane_distance - dot( ray_origin, plane_normal ) ) / denom;
 	return *t >= 0;
 }
 
 void barycentric( v3 a, v3 b, v3 c, v3 p, float * u, float * v, float * w ) {
 	v3 ab = b - a;
 	v3 ac = c - a;
 	v3 ap = p - a;
 
 	float dabab = dot( ab, ab );
 	float dabac = dot( ab, ac );
 	float dacac = dot( ac, ac );
 	float dapab = dot( ap, ab );
 	float dapac = dot( ap, ac );
 	float denom = dabab * dacac - dabac * dabac;
 
 	*v = ( dacac * dapab - dabac * dapac ) / denom;
 	*w = ( dabab * dapac - dabac * dapab ) / denom;
 	*u = 1.0f - *v - *w;
 }
 
 bool point_in_triangle( v3 p, v3 a, v3 b, v3 c ) {
 	float u, v, w;
 	barycentric( a, b, c, p, &u, &v, &w );
 
 	const float epsilon = 0.0001f;
 	if( u < -epsilon || u > 1.0f + epsilon ) return false;
 	if( v < -epsilon || v > 1.0f + epsilon ) return false;
 	if( w < -epsilon || w > 1.0f + epsilon ) return false;
 
 	return true;
 }
 
 // TODO: not robust w/ thin triangles
 v3 triangle_normal( v3 a, v3 b, v3 c ) {
 	return normalize( cross( b - a, c - a ) );
 }
 
 bool ray_vs_triangle( v3 ray_origin, v3 ray_dir, v3 p0, v3 p1, v3 p2, float * t ) {
 	v3 e01 = p1 - p0;
 	v3 e02 = p2 - p0;
 	v3 p = cross( ray_dir, e02 );
 	float det = dot( e01, p );
 
 	const float epsilon = 0.001f;
 	if( fabsf( det ) < epsilon ) return false;
 
 	float inv_det = 1.0f / det;
 
 	v3 d = ray_origin - p0;
 	float u = dot( d, p ) * inv_det;
 	if( u < 0 || u > 1 ) return false;
 
 	v3 q = cross( d, e01 );
 	float v = dot( ray_dir, q ) * inv_det;
 	if( v < 0 || u + v > 1 ) return false;
 
 	*t = dot( e02, q ) * inv_det;
 	return *t >= 0;
 }
 
 static size_t child_idx( size_t node_idx, size_t quadrant ) {
 	return node_idx * 4 + quadrant + 1;
 }
 
 static void sort4( int * indices, float * elems ) {
 	for( int i = 0; i < 4; i++ )
 		indices[ i ] = i;
 
 	// sort a/b and c/d
 	if( elems[ indices[ 0 ] ] > elems[ indices[ 1 ] ] )
 		swap2( indices[ 0 ], indices[ 1 ] );
 	if( elems[ indices[ 2 ] ] > elems[ indices[ 3 ] ] )
 		swap2( indices[ 2 ], indices[ 3 ] );
 
 	// find lowest and highest
 	if( elems[ indices[ 0 ] ] > elems[ indices[ 2 ] ] )
 		swap2( indices[ 0 ], indices[ 2 ] );
 	if( elems[ indices[ 1 ] ] > elems[ indices[ 3 ] ] )
 		swap2( indices[ 1 ], indices[ 3 ] );
 
 	// sort middle two
 	if( elems[ indices[ 1 ] ] > elems[ indices[ 2 ] ] )
 		swap2( indices[ 1 ], indices[ 2 ] );
 }
 
 static bool ray_vs_quadtree_node(
 	const QuadTree & qt, size_t node_idx, MinMaxu32 aabb,
 	const Ray3 & ray, float * t, v3 * xnormal, float max_dist
 ) {
 	if( aabb.maxs.x - aabb.mins.x == 2 || aabb.maxs.y - aabb.mins.y == 2 ) {
 		bool hit = false;
 
 		for( size_t i = 0; i < 2; i++ ) {
 			for( size_t j = 0; j < 2; j++ ) {
 				u32 x = aabb.mins.x + j;
 				u32 y = aabb.mins.y + i;
 
 				v3 p0 = heightmap_point( qt.heightmap, x + 0, y + 0 );
 				v3 p1 = heightmap_point( qt.heightmap, x + 1, y + 0 );
 				v3 p2 = heightmap_point( qt.heightmap, x + 0, y + 1 );
 				v3 p3 = heightmap_point( qt.heightmap, x + 1, y + 1 );
 
 				// bottom right triangle
 				float t0;
 				v3 normal0 = triangle_normal( p0, p1, p3 );
 				bool hit0 = ray_vs_triangle( ray.origin, ray.dir, p0, p1, p3, &t0 );
 				if( hit0 ) {
 					if( !hit || t0 < *t ) {
 						*t = t0;
 						*xnormal = normal0;
 						hit = true;
 					}
 				}
 
 				// top left triangle
 				float t1;
 				v3 normal1 = triangle_normal( p0, p3, p2 );
 				bool hit1 = ray_vs_triangle( ray.origin, ray.dir, p0, p3, p2, &t1 );
 				if( hit1 ) {
 					if( !hit || t1 < *t ) {
 						*t = t1;
 						*xnormal = normal1;
 						hit = true;
 					}
 				}
 			}
 		}
 
 		return hit;
 	}
 
 	MinMaxu32 aabbs[ 4 ];
 	float ts[ 4 ];
 	for( int i = 0; i < 4; i++ ) {
 		QuadTreeNode child = qt.nodes[ child_idx( node_idx, i ) ];
 		aabbs[ i ] = aabb.quadrant( i ).clamp_z( child.min_z, child.max_z );
 		ts[ i ] = FLT_MAX;
 		ray_vs_aabb( MinMax( aabbs[ i ] ), ray, &ts[ i ] );
 	}
 
 	int sorted[ 4 ];
 	sort4( sorted, ts );
 
 	for( int i = 0; i < 4; i++ ) {
 		int idx = sorted[ i ];
 		if( ts[ idx ] >= max_dist )
 			break;
 
 		if( ray_vs_quadtree_node( qt, child_idx( node_idx, idx ), aabbs[ idx ], ray, t, xnormal, max_dist ) ) {
 			return true;
 		}
 	}
 
 	return false;
 }
 
 bool ray_vs_quadtree( const QuadTree & qt, const Ray3 & ray, float * t, v3 * xnormal, float max_dist ) {
 	v3u32 mins = v3u32( 0, 0, qt.nodes[ 0 ].min_z );
 	v3u32 maxs = v3u32( qt.heightmap.w, qt.heightmap.h, qt.nodes[ 0 ].max_z );
 	MinMaxu32 aabb( mins, maxs );
 
 	return ray_vs_quadtree_node( qt, 0, aabb, ray, t, xnormal, max_dist );
 }
 
 bool ray_vs_terrain( const QuadTree & qt, const Ray3 & ray, float * t, v3 * xnormal, float max_dist ) {
 	float dont_care_t;
 	v3 dont_care_normal;
 
 	if( t == NULL )
 		t = &dont_care_t;
 	if( xnormal == NULL )
 		xnormal = &dont_care_normal;
 
 	bool hit_quadtree = ray_vs_quadtree( qt, ray, t, xnormal, max_dist );
 	if( hit_quadtree )
 		return true;
 
 	*xnormal = v3( 0, 0, 1 );
 	return ray_vs_plane( ray.origin, ray.dir, *xnormal, 0, t );
 }
 
 bool segment_vs_terrain( const QuadTree & qt, const v3 & start, const v3 & end, float * t, v3 * xnormal ) {
 	if( t != NULL )
 		*t = 1.0f;
 
 	if( start == end )
 		return false;
 
 	Ray3 ray( start, normalize( end - start ) );
 	float l = length( end - start );
 
 	float absolute_t;
 	bool ok = ray_vs_terrain( qt, ray, &absolute_t, xnormal, l );
 	if( !ok || absolute_t >= l )
 		return false;
 
 	if( t != NULL )
 		*t = absolute_t / l;
 	return true;
 }
 
 void bc5_to_heightmap( MemoryArena * arena, array2d< u16 > hm, u8 * bc5 ) {
 	MEMARENA_SCOPED_CHECKPOINT( arena );
 	array2d< v2 > rg = alloc_array2d< v2 >( arena, hm.w, hm.h );
 	decompress_bc5( rg, bc5 );
 
 	for( u32 y = 0; y < rg.h; y++ ) {
 		for( u32 x = 0; x < rg.w; x++ ) {
 			float r = rg( x, y ).x;
 			float g = rg( x, y ).y;
 			hm( x, y ) = u16( ( r * 256.0f + g ) * 255.0f );
 		}
 	}
 }