heightmap.c

#include <raylib.h>
#include <raymath.h>
#include <stdint.h>

#define min(a, b) (a) < (b) ? (a) : (b)
#define max(a, b) (a) > (b) ? (a) : (b)

static void
GridFromHeightMap(Vector3 *grid_positions, // output grid
                  int32_t  grid_width,     // point count x axis
                  int32_t  grid_length,    // point count z axis
                  float    scaleX,         // quad scale x
                  float    scaleZ,         // quad scale z
                  Color   *height_map,
                  int32_t  hm_width,
                  int32_t  hm_height,
                  float    hm_strength)
{
    float factor_x = (float)hm_width / grid_width;
    float factor_z = (float)hm_height / grid_length;

    for (int32_t gpz = 0; gpz < grid_length; ++gpz) {
        float hmy = floorf(factor_z * gpz);
        for (int32_t gpx = 0; gpx < grid_width; ++gpx) {
            float   hmx    = floorf(factor_x * gpx);
            float   height = height_map[(int32_t)(hmy * hm_width + hmx)].r / 255.0f * 2.0f - 1.0f;
            int32_t gp_idx = gpz * (grid_width) + gpx;

            grid_positions[gp_idx].x = scaleX * gpx;
            grid_positions[gp_idx].y = height * hm_strength;
            grid_positions[gp_idx].z = scaleZ * gpz;
        }
    }
}

// NOTE: Calculates 1 normal per point, averages surrounding "quad" normals for smooth shading
static void
CalculateGridNormals(const Vector3 *positions, Vector3 *normals, int32_t res_x, int32_t res_y)
{
    const Vector3 *v     = positions;
    int32_t        count = res_x * res_y;

    for (int32_t y = 0; y < res_y; ++y) {
        for (int32_t x = 0; x < res_x; ++x) {

            int32_t idx = y * (res_x) + x;

            Vector3 center = v[idx];
            Vector3 up     = v[max(0, y - 1) * (res_x) + x];
            Vector3 left   = v[max(0, idx - 1)];
            Vector3 right  = v[min(count - 1, idx + 1)];
            Vector3 down   = v[min(res_y - 1, y + 1) * (res_x) + x];

            Vector3 tan0 = Vector3Subtract(up, center);
            Vector3 bi0  = Vector3Subtract(left, center);
            Vector3 tan1 = Vector3Subtract(down, center);
            Vector3 bi1  = Vector3Subtract(right, center);

            Vector3 normal = Vector3CrossProduct(tan0, bi0);
            Vector3 n_mag  = Vector3Add(n_mag, normal);
          
            normal = Vector3CrossProduct(tan1, bi1);
            n_mag  = Vector3Add(n_mag, normal);

            normal = Vector3CrossProduct(bi1, tan0);
            n_mag  = Vector3Add(n_mag, normal);

            normal = Vector3CrossProduct(bi0, tan1);
            n_mag   = Vector3Add(n_mag, normal);

            normals[idx] = Vector3Normalize(n_mag);
        }
    }
}

// NOTE: Calculates normals based on mesh triangles. Does no smoothing
static void
CalculateVertexNormals(const Vector3 *v, Vector3 *n, int32_t count)
{
    for (int32_t i = 0; i < count; i += 3) {
        Vector3 tan    = Vector3Subtract(v[i + 1], v[i]);
        Vector3 bi     = Vector3Subtract(v[i + 2], v[i]);
        Vector3 normal = Vector3Normalize(Vector3CrossProduct(tan, bi));

        n[i + 0] = normal;
        n[i + 1] = normal;
        n[i + 2] = normal;
    }
}

static void
CalculatePlaneFromGrid(const Vector3 *grid_positions,
                       const Vector3 *grid_normals,
                       Vector3       *vertices,
                       Vector3       *normals,
                       int32_t        res_x,
                       int32_t        res_y)
{
    Vector3 *v = vertices;
    Vector3 *n = normals;

    for (int32_t y = 0; y < res_y; ++y) {
        for (int32_t x = 0; x < res_x; ++x) {
            // Read grid positions as "Quads". So four points for 2 triangles. 6 Tri vertices per loop.
            // NOTE: the grid's X and Y res is +1 so we are not reading off the end

            Vector3 v0 = grid_positions[(y + 0) * (res_x + 1) + (x + 1)];
            Vector3 v1 = grid_positions[(y + 0) * (res_x + 1) + (x + 0)];
            Vector3 v2 = grid_positions[(y + 1) * (res_x + 1) + (x + 0)];
            Vector3 v3 = grid_positions[(y + 1) * (res_x + 1) + (x + 1)];

            Vector3 n0 = grid_normals[(y + 0) * (res_x + 1) + (x + 1)];
            Vector3 n1 = grid_normals[(y + 0) * (res_x + 1) + (x + 0)];
            Vector3 n2 = grid_normals[(y + 1) * (res_x + 1) + (x + 0)];
            Vector3 n3 = grid_normals[(y + 1) * (res_x + 1) + (x + 1)];

            // Tri 1
            v[0] = v0;
            v[1] = v1;
            v[2] = v2;

            n[0] = n0;
            n[1] = n1;
            n[2] = n2;

            v += 3;
            n += 3;

            // Tri 2
            v[0] = v0;
            v[1] = v2;
            v[2] = v3;

            n[0] = n0;
            n[1] = n2;
            n[2] = n3;

            v += 3;
            n += 3;
        }
    }
}

static Mesh
CreatePlaneMesh(int32_t res_x, int32_t res_y)
{
    Mesh mesh          = {0};
    mesh.triangleCount = 2 * res_x * res_y;
    mesh.vertexCount   = mesh.triangleCount * 3;
    mesh.vertices      = MemAlloc(mesh.vertexCount * sizeof(float) * 3);
    mesh.texcoords     = MemAlloc(mesh.vertexCount * sizeof(float) * 2);
    mesh.normals       = MemAlloc(mesh.vertexCount * sizeof(float) * 3);

    // Calculate UVs
    Vector2 *tc    = (Vector2 *)mesh.texcoords;
    int32_t  count = mesh.vertexCount;
    for (int32_t y = 0; y < res_y; ++y) {
        for (int32_t x = 0; x < res_x; ++x) {
            float rx = 1.0f / (float)res_x;
            float ry = 1.0f / (float)res_y;

            float uvx_min = rx * x;
            float uvy_min = ry * y;

            float uvx_max = uvx_min + rx;
            float uvy_max = uvy_min + ry;

            tc[0].x = uvx_max, tc[0].y = uvy_min;
            tc[1].x = uvx_min, tc[1].y = uvy_min;
            tc[2].x = uvx_min, tc[2].y = uvy_max;

            tc += 3;

            tc[0].x = uvx_max, tc[0].y = uvy_min;
            tc[1].x = uvx_min, tc[1].y = uvy_max;
            tc[2].x = uvx_max, tc[2].y = uvy_max;

            tc += 3;
        }
    }

    UploadMesh(&mesh, false);
    return mesh;
}

int
main(void)
{
    SetConfigFlags(FLAG_MSAA_4X_HINT);
    InitWindow(1920, 1080, "Heightmap");
    SetTargetFPS(60);
    DisableCursor();

    Image  height_map        = GenImagePerlinNoise(512, 512, 0, 0, 3.0f);
    Color *height_map_colors = LoadImageColors(height_map);

    int32_t plane_res_x = 64;
    int32_t plane_res_y = 64;
    Mesh    plane       = CreatePlaneMesh(plane_res_x, plane_res_y);

    // NOTE: +1 to each dimension because we are making QUADS so we need +1 to form the last row/column
    int32_t  grid_point_count_x = plane_res_x + 1;
    int32_t  grid_point_count_y = plane_res_y + 1;
    int32_t  grid_point_count   = grid_point_count_x * grid_point_count_y;
    Vector3 *grid_positions     = MemAlloc(grid_point_count * sizeof(Vector3));
    Vector3 *grid_normals       = MemAlloc(grid_point_count * sizeof(Vector3));

    float width_meters  = 10.0f / plane_res_x;
    float length_meters = 10.0f / plane_res_y;

    float height_map_strength = 0.5f;

    GridFromHeightMap(grid_positions,
                      grid_point_count_x,
                      grid_point_count_y,
                      width_meters,
                      length_meters,
                      height_map_colors,
                      height_map.width,
                      height_map.height,
                      height_map_strength);

    CalculateGridNormals(grid_positions, grid_normals, grid_point_count_x, grid_point_count_y);
    CalculatePlaneFromGrid(grid_positions, grid_normals, (Vector3 *)plane.vertices, (Vector3 *)plane.normals, plane_res_x, plane_res_y);

    // calculates non-smooth normals for each triangle
    // CalculateVertexNormals((Vector3 *)plane.vertices, (Vector3 *)plane.normals, plane.vertexCount);

    Model model = LoadModelFromMesh(plane);

    UpdateMeshBuffer(plane, 0, plane.vertices, sizeof(Vector3) * plane.vertexCount, 0);
    UpdateMeshBuffer(plane, 1, plane.texcoords, sizeof(Vector2) * plane.vertexCount, 0);
    UpdateMeshBuffer(plane, 2, plane.normals, sizeof(Vector3) * plane.vertexCount, 0);

    Camera camera     = {0};
    camera.position   = (Vector3){0.0f, 3.0f, 5.0f};
    camera.target     = (Vector3){0.3f, 0.0f, 3.0f};
    camera.up         = (Vector3){0.0f, 1.0f, 0.0f};
    camera.fovy       = 45.0f;
    camera.projection = CAMERA_PERSPECTIVE;

    while (!WindowShouldClose()) {

        UpdateCamera(&camera, CAMERA_FREE);

        BeginDrawing();
        ClearBackground(RAYWHITE);
        BeginMode3D(camera);

        // DrawSphereWires(Vector3Zero(), 2.0f, 6, 6, YELLOW);
        DrawModelWires(model, Vector3Zero(), 1.0f, GREEN);
        DrawModel(model, Vector3Zero(), 1.0f, DARKGRAY);

        EndMode3D();

        EndDrawing();
    }

    MemFree(grid_positions);
    MemFree(grid_normals);

    UnloadModel(model);

    UnloadImageColors(height_map_colors);
    UnloadImage(height_map);
    CloseWindow();
}