projecting 3d points to 2d in webassembly

proj.c

// Copyright (c) 2021, Wei Mingzhi.
// Released under BSD License.
// https://opensource.org/licenses/BSD-3-Clause

/*
 compile with:

 clang --target=wasm32 \
  -O3 \
    -flto -nostdlib \
    -Wl,--no-entry -Wl,--export-all -Wl,--lto-O3 \
    -Wl,-z,stack-size=512000 -Wl,--import-memory -o proj.wasm proj.c
*/

extern unsigned char __heap_base;

unsigned char *bump_pointer = &__heap_base;

void* allocMem(int n) {
    void *r = bump_pointer;
    bump_pointer += n;
    return r;
}

void resetMem() {
    bump_pointer = &__heap_base;
}

typedef float vec3_t[3];
typedef float mat4_t[16];

static inline void vec_applyMatrix4(vec3_t v, mat4_t e) {
    const float x = v[0], y = v[1], z = v[2];
    const float w = 1 / (e[3] * x + e[7] * y + e[11] * z + e[15]);

    v[0] = (e[0] * x + e[4] * y + e[8] * z + e[12]) * w;
    v[1] = (e[1] * x + e[5] * y + e[9] * z + e[13]) * w;
    v[2] = (e[2] * x + e[6] * y + e[10] * z + e[14]) * w;
}

static inline void mat_multiplyMatrices(mat4_t ae, mat4_t be, mat4_t te) {
    const float a11 = ae[0], a12 = ae[4], a13 = ae[8], a14 = ae[12];
    const float a21 = ae[1], a22 = ae[5], a23 = ae[9], a24 = ae[13];
    const float a31 = ae[2], a32 = ae[6], a33 = ae[10], a34 = ae[14];
    const float a41 = ae[3], a42 = ae[7], a43 = ae[11], a44 = ae[15];
    const float b11 = be[0], b12 = be[4], b13 = be[8], b14 = be[12];
    const float b21 = be[1], b22 = be[5], b23 = be[9], b24 = be[13];
    const float b31 = be[2], b32 = be[6], b33 = be[10], b34 = be[14];
    const float b41 = be[3], b42 = be[7], b43 = be[11], b44 = be[15];

    te[0] = a11 * b11 + a12 * b21 + a13 * b31 + a14 * b41;
    te[4] = a11 * b12 + a12 * b22 + a13 * b32 + a14 * b42;
    te[8] = a11 * b13 + a12 * b23 + a13 * b33 + a14 * b43;
    te[12] = a11 * b14 + a12 * b24 + a13 * b34 + a14 * b44;

    te[1] = a21 * b11 + a22 * b21 + a23 * b31 + a24 * b41;
    te[5] = a21 * b12 + a22 * b22 + a23 * b32 + a24 * b42;
    te[9] = a21 * b13 + a22 * b23 + a23 * b33 + a24 * b43;
    te[13] = a21 * b14 + a22 * b24 + a23 * b34 + a24 * b44;

    te[2] = a31 * b11 + a32 * b21 + a33 * b31 + a34 * b41;
    te[6] = a31 * b12 + a32 * b22 + a33 * b32 + a34 * b42;
    te[10] = a31 * b13 + a32 * b23 + a33 * b33 + a34 * b43;
    te[14] = a31 * b14 + a32 * b24 + a33 * b34 + a34 * b44;

    te[3] = a41 * b11 + a42 * b21 + a43 * b31 + a44 * b41;
    te[7] = a41 * b12 + a42 * b22 + a43 * b32 + a44 * b42;
    te[11] = a41 * b13 + a42 * b23 + a43 * b33 + a44 * b43;
    te[15] = a41 * b14 + a42 * b24 + a43 * b34 + a44 * b44;
}

typedef struct {
    vec3_t normal;
    float constant;
} plane_t;

static inline void plane_setComponents(plane_t *plane, float x, float y, float z, float w) {
    plane->normal[0] = x;
    plane->normal[1] = y;
    plane->normal[2] = z;
    plane->constant = w;
}

#define FLT_MAX 3.402823466e+38F

static inline float Q_rsqrt(float number) {
    if (number <= 1e-10f) {
        return FLT_MAX;
    }
    return 1.0f / sqrtf(number);
}

static inline void plane_normalize(plane_t *plane) {
    float inverseNormalLength = Q_rsqrt(plane->normal[0] * plane->normal[0] +
                                        plane->normal[1] * plane->normal[1] +
                                        plane->normal[2] * plane->normal[2]);
    plane->normal[0] *= inverseNormalLength;
    plane->normal[1] *= inverseNormalLength;
    plane->normal[2] *= inverseNormalLength;
    plane->constant *= inverseNormalLength;
}

static inline float plane_distanceToPoint(plane_t *plane, vec3_t point) {
    return plane->normal[0] * point[0] +
           plane->normal[1] * point[1] +
           plane->normal[2] * point[2] +
           plane->constant;
}

typedef struct {
    plane_t planes[6];
} frustum_t;

static inline void frustum_setFromProjectionMatrix(frustum_t *f, mat4_t me) {
    float me0 = me[0], me1 = me[1], me2 = me[2], me3 = me[3];
    float me4 = me[4], me5 = me[5], me6 = me[6], me7 = me[7];
    float me8 = me[8], me9 = me[9], me10 = me[10], me11 = me[11];
    float me12 = me[12], me13 = me[13], me14 = me[14], me15 = me[15];

    plane_setComponents(&(f->planes[0]), me3 - me0, me7 - me4, me11 - me8, me15 - me12);
    plane_setComponents(&(f->planes[1]), me3 + me0, me7 + me4, me11 + me8, me15 + me12);
    plane_setComponents(&(f->planes[2]), me3 + me1, me7 + me5, me11 + me9, me15 + me13);
    plane_setComponents(&(f->planes[3]), me3 - me1, me7 - me5, me11 - me9, me15 - me13);
    plane_setComponents(&(f->planes[4]), me3 - me2, me7 - me6, me11 - me10, me15 - me14);
    plane_setComponents(&(f->planes[5]), me3 + me2, me7 + me6, me11 + me10, me15 + me14);
    plane_normalize(&(f->planes[0]));
    plane_normalize(&(f->planes[1]));
    plane_normalize(&(f->planes[2]));
    plane_normalize(&(f->planes[3]));
    plane_normalize(&(f->planes[4]));
    plane_normalize(&(f->planes[5]));
}

static inline int frustum_containsPoint(frustum_t *f, vec3_t point) {
    for (int i = 0; i < 6; i++) {
        float dist = plane_distanceToPoint(&(f->planes[i]), point);
        if (dist < 0) {
            return 0;
        }
    }
    return 1;
}

static frustum_t frustum;

void doProject(vec3_t points, int pointNum, mat4_t matrixWorldInverse, mat4_t projectionMatrix, float w, float h) {
    mat4_t frustumMat;
    mat_multiplyMatrices(projectionMatrix, matrixWorldInverse, frustumMat);
    frustum_setFromProjectionMatrix(&frustum, frustumMat);
    h *= -1;

    for (int i = 0; i < pointNum; i++) {
        float *vec = &points[i * 3];
        if (frustum_containsPoint(&frustum, vec)) {
            vec_applyMatrix4(vec, frustumMat);
            vec[0] = (int)(vec[0] * w + w + 0.5f);
            vec[1] = (int)(vec[1] * h - h + 0.5f);
            vec[2] = 1;
        } else {
            vec[2] = 0;
        }
    }
}