Terrain interpolation algorithm

can unity release your terrain vertex height interpolation algorithm? i need to do player’s position validation on server side which need the acurate height where the player stand. I can export the heightmap from unity, but it’s not a per vertex heightmap, so interpolation is needed.

Or can some one give me some other way to reach the goal?

You could do raycast down?

Or this,

My application server is not a unity3d project(it should be, right?), i can’t call unity’s api at all. i can only import the unity3d heightmap and do the job by my own code.

I’d like to see an optimized algorithm also. I have a similar situation where my C# 4.0 server has the terrain height data, and has control over player and item positioning. The server needs to have the final say over where something rests on the terrain, without relying on Unity.

Any one can give out a solution?

just use bi-linear interpolation algorithm its basic but works fine. basically we are interpolating the unit square with respect to x and z(boo… calculus :p)

var x = Random.Range(0, Terrain.activeTerrain.terrainData.size.x);
var z = Random.Range(0, Terrain.activeTerrain.terrainData.size.z);

var hx0z0 = Terrain.activeTerrain.SampleHeight(new Vector3(Mathf.FloorToInt(x), 0, Mathf.FloorToInt(z))); // known height (x0, z0)
var hx1z0 = Terrain.activeTerrain.SampleHeight(new Vector3(Mathf.CeilToInt(x), 0, Mathf.FloorToInt(z))); // known height (x1, z0)
var hx0z1 = Terrain.activeTerrain.SampleHeight(new Vector3(Mathf.FloorToInt(x), 0, Mathf.CeilToInt(z))); // known height (x0, z1)
var hx1z1 = Terrain.activeTerrain.SampleHeight(new Vector3(Mathf.CeilToInt(x), 0, Mathf.CeilToInt(z))); // known height (x1, z1)

var u0v0 = hx0z0 * (Mathf.CeilToInt(x) - x) * (Mathf.CeilToInt(z) - z); // interpolated (x0, z0)
var u1v0 = hx1z0 * (x - Mathf.FloorToInt(x)) * (Mathf.CeilToInt(z) - z); // interpolated (x1, z0)
var u0v1 = hx0z1 * (Mathf.CeilToInt(x) - x) * (z - Mathf.FloorToInt(z)); // interpolated (x0, z1)
var u1v1 = hx1z1 * (x - Mathf.FloorToInt(x)) * (z - Mathf.FloorToInt(z)); // interpolated (x1, z1)

var h = u0v0 + u1v0 + u0v1 + u1v1; // estimate

Or you can use the fast non-linear algorithm which is basically the product of x and z interpolation.

var h = hx0z0 + (hx1z0 - hx0z0) * (x - (int)x) + (hx0z1 - hx0z0) * (z - (int)z) + (hx0z0 - hx1z0 - hx0z1 + hx1z1) * (x - (int)x) * (z - (int)z);

I used Terrain.SampleHeight(x, z) to get the known height at fixed positions but in your server you might have an array of terrain heights, so you can use them to get the height at known positions then use the same code to interpolate between them

1 Like

Thank you that was very useful, here is a version for use with texture files, I found a bizarre issue with it, if I use a simple divide for the height map for example X/10, Y/10, to make the map 10 times wider, every time it finds a vertex position X and Y which has no decimal numbers, in other words and then to get, it gets it wrong and returns the result is 0. That is why at the top of this code there is a funny multiplier for X and Y map texture positions.

function  equation1  ( vertex: Vector3 ): float {// texture reading-w.interpolation

    var  x = vertex.x*.30101+.1;
    var  z = vertex.z*.30101+.1;


 

var hx0z0 = mytexture.GetPixel(Mathf.FloorToInt(x), Mathf.FloorToInt(z)).grayscale; // known height (x0, z0)

var hx1z0 = mytexture.GetPixel(Mathf.CeilToInt(x), Mathf.FloorToInt(z)).grayscale; // known height (x1, z0)

var hx0z1 = mytexture.GetPixel(Mathf.FloorToInt(x), Mathf.CeilToInt(z)).grayscale; // known height (x0, z1)

var hx1z1 = mytexture.GetPixel(Mathf.CeilToInt(x), Mathf.CeilToInt(z)).grayscale; // known height (x1, z1)

 

var u0v0 = hx0z0 * (Mathf.CeilToInt(x) - x) * (Mathf.CeilToInt(z) - z); // interpolated (x0, z0)

var u1v0 = hx1z0 * (x - Mathf.FloorToInt(x)) * (Mathf.CeilToInt(z) - z); // interpolated (x1, z0)

var u0v1 = hx0z1 * (Mathf.CeilToInt(x) - x) * (z - Mathf.FloorToInt(z)); // interpolated (x0, z1)

var u1v1 = hx1z1 * (x - Mathf.FloorToInt(x)) * (z - Mathf.FloorToInt(z)); // interpolated (x1, z1)

 

var h = u0v0 + u1v0 + u0v1 + u1v1; // estimate	
//var	h =  mytexture.GetPixel(vertex.x*.51, vertex.z*.51).grayscale * 88;
  //  if ( vertex.z == 0 ){ Debug.Log(  " pseudo-  "+vertex.x +  " result "   +vertex.y);}		
    return h*88;

}

<a href="http://tinypic.com?ref=5czqzl" target="_blank"><img src="http://i48.tinypic.com/5czqzl.jpg" border="0" alt="Image and video hosting by TinyPic"></a>There is one small problem though-on slopes that are slanted in both X and Y, the interpolation seems to produce errors in between the pixels, it looks like mottled instead of straight, I don’t think that is due to the source texture, I think the interpolation is making some mottled results. can you explain that?

So, while it has been many years, the problem remains and is getting more relevant by a day. Especially now when DOTS coming up and there is no way to estimate interpolated heights manually. It is also noticeable that this thread is coming to one of the top results when searching for Unity heightmap interpolation.

But this probably will be no more as a problem. I was recently looking a bit more into details and made some comparison tests here GitHub - chanfort/Terrain-heightmap-manual-interpolation.

It looks like the errors ZoomDomain got might be explained that terrain heightmap is not just a grid of quads. Instead, each quad is represented as two triangles with the diagonal going through the middle in (0,0) => (1,1) direction. As a result, terrain can bend through these diagonals producing errors. As a result I looked a bit more into a different barycentric interpolation and added solution code (based on http://wiki.unity3d.com/index.php?title=Barycentric&oldid=19264) in the repo. It seems like maximum error over large set of points is decreasing close to floating point precision limit which may mean that it is the correct approach to use.

I also added jobified and bursted code example and just a side note that it outperforms regular calculation by two orders of magnitude (almost 100x times).

Just a reminder that getting interpolated terrain heights is essential when placing objects on the surface of terrain. Applications of this ranges between spawning trees, grass, detail objects as well as adjusting heights for NavMeshAgents in order to walk on the terrain at the correct height (i.e. preventing walking in the air or having sinked legs where NavMesh does not match terrain).

3 Likes

Another approach I found that works really well is the code near the bottom of this page by @joshrs926 : How is Terrain.SampleHeight implemented? - Questions & Answers - Unity Discussions . It uses triangles to do the interpolation, but I’m not sure how fast or accurate it is compared to @chanfort 's barycentric method.

It was very easy to add to my code - I changed the struct to an IComponentData with a blob-array-ref and then I can simply call dstManager.AddComponentData(entity, new TerrainHeightData(terrain)); (where the terrain parameter is a UnityEngine.Terrain). Then my terrain entity will have a fancy burst-and-thread-safe terrain data component that I can use to fetch terrain heights like this:

TerrainHeightData terrainMap = GetSingleton<TerrainHeightData>();
float height = terrainMap.SampleHeight(myPos);

If anyone wants my blobarray version of the code, just let me know and I’ll share it - 90% of it corresponds with the code in the thread I posted above.

For things that don’t happen very often, like the player placing objects on a terrain, raycasting works just fine and is trivial to implement. But if you have thousands of navmesh agents or something and they all need to fetch an altitude every frame, then it’s probably smart to avoid raycasting and fetch the altitude directly via one of the methods in this thread.

I’m glad it worked for you! I bet @chanfort method is good too since he recognized that Unity appears to represent a terrain as a grid of quads each being split into 2 tris with their border going from bottom left to top right (at least at the lowest terrain LOD). I bet his or my method are both accurate though I’m not sure whose is faster. As a side note, I was actually later able to replace the array of floats with an array of ushorts to halve the memory requirement. I also implemented holes as well. You can pass in a UnityEngine.Terrain and my DOTS terrain copy will be able to tell you if a certain sample position is in a hole or not. Perhaps I’ll post that code on that thread at some point.

1 Like

Converting to ushort is a great idea. I just did it and I can’t tell the difference. I’m not using holes at the moment, but maybe later I’ll need it, so if you post it, I hope you do it in this thread so I can get notified :-). Thanks so much for sharing that script - it saved me a ton of headache!

One last optimization… There is a division operation when converting from ushort and there are two more division operations that get called 3x per lookup (so 7 total) which can be converted to the shift operation to improve the speed even more. There are still two divisions per sample-lookup that can’t be converted, but overall 2 is a lot better than 9.

That’s interesting! I’ll make you a deal. If you post your shift operation thing then I’ll post my holes thing! Have you tested any performance increase with that shift operation?

I haven’t tested the performance of this particular class, but I know from reading and prior testing that shifting is much faster than divide. I’m using the terrain sampling to get the height above terrain for tens of thousands of agents once per frame, so I suspect that eliminating as many divides as possible is helpful in that scenario.

// Most of this code came from joshrs926 at https://answers.unity.com/questions/1642900/how-is-terrainsampleheight-implemented.html
// I modified it to be an IComponentData and converted the height map from an unsafe-list to a blob-array.
// I made it store terrain as 16 bit ushorts instead of 32 bit floats to save memory.
// I changed some division operations to shifts to make sample operations run faster.
//
// To place it on a terrain during a conversion, use:
//    Terrain terrain = GetComponent<Terrain>();
//    dstManager.AddComponentData(entity, new TerrainHeightData(terrain));
//
// To read height data in a job:
//    TerrainHeightData terrainMap = GetSingleton<TerrainHeightData>()
//    float height = terrainMap.SampleHeight(pos);

using Unity.Collections;
using Unity.Mathematics;
using UnityEngine;
using Unity.Entities;
    
public struct TerrainHeightData : IComponentData
{
    public AABB AABB { get; private set; }
    public bool IsValid => heightMapRef.IsCreated;

    public BlobAssetReference<BlobArray<ushort>> heightMapRef;

    int resolution;
    float2 sampleSize;
    int QuadCount => resolution - 1;

    float minVal; // used to convert to/from 16 bit
    float maxVal; // used to convert to/from 16 bit
    
    // Constructor that fills all the variables (including height-map blob ref), given the specified terrain object.
    public TerrainHeightData(Terrain terrain)
    {
        resolution = terrain.terrainData.heightmapResolution;
        sampleSize = new float2(terrain.terrainData.heightmapScale.x, terrain.terrainData.heightmapScale.z);
        AABB = GetTerrrainAABB(terrain);

        using BlobBuilder builder = new BlobBuilder(Allocator.Temp);
        ref BlobArray<ushort> root = ref builder.ConstructRoot<BlobArray<ushort>>();
        BlobBuilderArray<ushort> arrayBuilder = builder.Allocate(ref root, resolution * resolution);
       
        var map = terrain.terrainData.GetHeights(0, 0, resolution, resolution);

        // First the min an max values need to be calculated and stored (for later use when decoding a height).
        // This allows the height map to be of type ushort, which is half the size of float.
        minVal = float.MaxValue;
        maxVal = float.MinValue;
        for (int y = 0; y < resolution; y++) {
            for (int x = 0; x < resolution; x++) {
                if (map[y, x] < minVal) { minVal = map[y,x]; }
                if (map[y, x] > maxVal) { maxVal = map[y,x]; }
            }
        }

        // now store each point, scaling it via the min/max calculated above in the process
        for (int y = 0; y < resolution; y++) {
            for (int x = 0; x < resolution; x++) {
                int i = y * resolution + x;
                arrayBuilder[i] = (ushort)MathUtils.Scale(map[y, x], minVal, maxVal, 0, ushort.MaxValue);
            }
        }
        heightMapRef = builder.CreateBlobAssetReference<BlobArray<ushort>>(Allocator.Persistent);
    }

    // Returns world height of terrain at x and z position values.
    public float SampleHeight(float3 worldPosition)
    {
        GetTriAtPosition(worldPosition, out Triangle tri);
        return tri.SampleHeight(worldPosition);
    }
    
    // Returns world height of terrain at x and z position values. Also outputs normalized normal vector of terrain at position.
    public float SampleHeight(float3 worldPosition, out float3 normal)
    {
        GetTriAtPosition(worldPosition, out Triangle tri);
        normal = tri.Normal;
        return tri.SampleHeight(worldPosition);
    }
    
    // fetches the triangle at the specified position
    void GetTriAtPosition(float3 worldPosition, out Triangle tri)
    {
        if (!IsWithinBounds(worldPosition)) {
            throw new System.ArgumentException($"Position given {worldPosition} is outside of terrain x or z bounds.");
        }
        float2 localPos = new float2(worldPosition.x - AABB.Min.x, worldPosition.z - AABB.Min.z);
        float2 samplePos = localPos / sampleSize;
        int2 sampleFloor = (int2)math.floor(samplePos);
        float2 sampleDecimal = samplePos - sampleFloor;
        bool upperLeftTri = sampleDecimal.y > sampleDecimal.x;
        int2 v1Offset = upperLeftTri ? new int2(0, 1) : new int2(1, 1);
        int2 v2Offset = upperLeftTri ? new int2(1, 1) : new int2(1, 0);
        float3 v0 = GetWorldVertex(sampleFloor);
        float3 v1 = GetWorldVertex(sampleFloor + v1Offset);
        float3 v2 = GetWorldVertex(sampleFloor + v2Offset);
        tri = new Triangle(v0, v1, v2);
    }
    
    //public void Dispose()
    //{
    //    heightMap.Dispose();
    //}
    
    // returns true if the specified point is within the terrain bounds (on the x/z plane)
    bool IsWithinBounds(float3 worldPos)
    {
        return worldPos.x >= AABB.Min.x && worldPos.z >= AABB.Min.z && worldPos.x <= AABB.Max.x && worldPos.z <= AABB.Max.z;
    }
    
    // return the vertex at the specified hight map coordinates.
    float3 GetWorldVertex(int2 heightMapCrds)
    {
        int i = heightMapCrds.x + heightMapCrds.y * resolution;
        // scale the height (which is a ushort) back to a float using the precalculated min/max vals
        // Since val is being divided by 65536, which is a power of 2, a fast divide can be used.
        float scaledHeight = MathUtils.FastDivideByPow2(heightMapRef.Value[i], 65536) * (maxVal - minVal) + minVal;
        float3 vertexPercentages = new float3(MathUtils.FastDivideByPow2(heightMapCrds.x, (uint)QuadCount), scaledHeight, MathUtils.FastDivideByPow2(heightMapCrds.y, (uint)QuadCount));
        return AABB.Min + AABB.Size * vertexPercentages;
    }
    
    // return the AABB box for the given terrain
    static AABB GetTerrrainAABB(Terrain terrain)
    {
        float3 min = terrain.transform.position;
        float3 max = min + (float3)terrain.terrainData.size;
        float3 extents = (max - min) * 0.5f;
        return new AABB() { Center = min + extents, Extents = extents };
    }
}

public readonly struct Triangle
{
    public float3 V0 { get; }
    public float3 V1 { get; }
    public float3 V2 { get; }
    // this is already normalized
    public float3 Normal { get; }
    
    public Triangle(float3 v0, float3 v1, float3 v2)
    {
        V0 = v0;
        V1 = v1;
        V2 = v2;
        Normal = math.normalize(math.cross(V1 - V0, V2 - V0));
    }
    public float SampleHeight(float3 position)
    {
        // plane formula: a(x - x0) + b(y - y0) + c(z - z0) = 0
        // <a,b,c> is a normal vector for the plane
        // (x,y,z) and (x0,y0,z0) are any points on the plane           
        return (-Normal.x * (position.x - V0.x) - Normal.z * (position.z - V0.z)) / Normal.y + V0.y;
    }
}
// this is collection of useful static math functions
//--------------------------------------------------------------------------------------------------//

using Unity.Mathematics;
using System.Runtime.CompilerServices;

public static class MathUtils
{
    // These mysterious kotay bits are used to access a bit mask used when calculating the power of 2 (log2(n)).
    static readonly byte[] KotayBits = { 0, 1, 2, 16, 3, 6, 17, 21, 14, 4, 7, 9, 18, 11, 22, 26, 31, 15, 5, 20, 13, 8, 10, 25, 30, 19, 12, 24, 29, 23, 28, 27 };

    // there are 31 slots in this array, where each slot corresponds with a power of 2.
    // the array contains the inverse of each power.
    static readonly float[] InversePowers = {
        1f, 1f/2f, 1f/4f, 1f/8f, 1f/16f, 1f/32f, 1f/64f, 1f/128f, 1f/256f, 1f/512f, 1f/1024f, 1f/2048f, 1f/4096f, 1f/8192f, 1f/16384f, 1f/32768f, 1f/65536f, 1f/131072f, 1f/262144f, 1f/524288f, 1f/1048576f,
        1f/2097152f, 1f/4194304f, 1f/8388608f, 1f/16777216f, 1f/33554432f, 1f/67108864f, 1f/134217728f, 1f/268435456f, 1f/536870912f, 1f/1073741824f, 1f/2147483648f, 1f/4294967296f
    };

    // WARNING - This ONLY works if the numerator is an unsigned int, and the divisor is a power of 2!.
    // WARNING - Also note that the return value is a uint - NOT a float or an int!.
    // See http://sree.kotay.com/2007/04/shift-registers-and-de-bruijn-sequences_10.html
    // It is slightly faster than the float version of this function (one less lookup and shift instead of multiply).
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static uint FastIntegerDivideByPow2(uint num, uint den)
    {
        //0x04ad19df --> 78453215 (special kotay bit mask).
        // the kotay bits give us a super fast way of doing a log2(n) (aka finding the power of 2 for n).
        return num >> KotayBits[den * 78453215U >> 27];
    }

    // WARNING - This ONLY works if the divisor is a power of 2!
    // den --> the largest possible value is 2^31 (2147483648).
    // The return value is a float.
    // It is slightly slower than the integer version of this function (extra lookup and multiply instead of shift).
    // See http://sree.kotay.com/2007/04/shift-registers-and-de-bruijn-sequences_10.html
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static float FastDivideByPow2(float num, uint den)
    {
        return InversePowers[KotayBits[den * 78453215U >> 27]] * num;
    }

    // same as above, except this accepts a float3
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static float3 FastDivideByPow2(float3 num, uint den)
    {
        return new float3(FastDivideByPow2(num.x, den), FastDivideByPow2(num.y, den), FastDivideByPow2(num.z, den));
    }

    // scales val, which is between oldMin and oldMax, to a number between newMin and newMax
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static float Scale(float val, float oldMin, float oldMax, float newMin, float newMax)
    {
        return (((val - oldMin) / (oldMax - oldMin)) * (newMax - newMin)) + newMin;
    }
}