Inverse distance weighting

idw.c

#include <assert.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

typedef struct {
	size_t size[2];
	size_t stride[2];
	float *elements;
} matrix_t;

float *mat_ij(const matrix_t *m, size_t row, size_t col) {
	assert(0 <= row && row < m->size[0]);
	assert(0 <= col && col < m->size[1]);
	return &(m->elements[row*m->stride[0] + col*m->stride[1]]);
}

size_t mat_rows(const matrix_t *m) { return m->size[0]; }
size_t mat_cols(const matrix_t *m) { return m->size[1]; }

void idw_layer(float *x, const matrix_t *nx, const matrix_t *ny, float *y) {
	assert(nx->size[0] == ny->size[0]);

	// clear y.
	for (size_t i = 0; i < mat_cols(ny); ++i) {
		y[i] = 0;
	}	

	float cum_w = 0;
	for (size_t r = 0; r < mat_rows(nx); ++r) {
		// compute distance to node.
		float dist = 0;
		for (size_t c = 0; c < mat_cols(nx); ++c) {
			float delta = x[c] - *mat_ij(nx, r,c);
			dist += delta*delta;
		}

		float w = 1.0f/fmaxf(1e-4, dist);

		// add unnormalized weighted to y.
		for (size_t c = 0; c < mat_cols(ny); ++c) {
			y[c] += w * *mat_ij(ny,r,c);
		}
		cum_w += w;
	}
	// normalize y.
	for (size_t i = 0; i < mat_cols(ny); ++i) {
		y[i] /= cum_w;
	}	
}

int main() {
	const int width = 100, height = 100;

	FILE *f = fopen("idw.pgm", "w");
	assert(f);


	matrix_t nx = {
		.size = {6,2},
		.stride = {2,1},
		.elements = (float[6*2]) {0}
	};
	matrix_t ny = {
		.size = {6,3},
		.stride = {3,1},
		.elements = (float[6*3]) {0}
	};

	// Randomly initialize IDW layer.
	srand(time(NULL));
	for (size_t i = 0; i < 6*2; ++i) {
		nx.elements[i] = rand() % 100;
	}
	for (size_t i = 0; i < 6*3; ++i) {
		ny.elements[i] = rand() % 255;
	}


	fprintf(f, "P3\n%d %d\n255\n", width, height);
	for (int x = 0; x < width; ++x) {
		for (int y = 0; y < height; ++y) {
			float pos[2] = {x, y};
			float rgb[3];
			idw_layer(pos, &nx, &ny, rgb);

			fprintf(f, "%d %d %d  ", (int) rgb[0], (int) rgb[1], (int) rgb[2]);
		}
		fprintf(f, "\n");
	}
}