Calculate edge-edge intersections of convex polygons via Sutherland-Hodgman in 2D

edge_interesections.cpp

#ifndef _CRT_SECURE_NO_WARNINGS
#	define _CRT_SECURE_NO_WARNINGS
#endif
#ifndef _CRT_NONSTDC_NO_DEPRECATE
#	define _CRT_NONSTDC_NO_DEPRECATE
#endif

#include <vector>

// -------------------------------------------------------------------------------------------------
// Basic 2D vector math.

struct v2
{
	v2() { }
	v2(float x, float y) { this->x = x; this->y = y; }
	float x;
	float y;
	v2 operator-() { return v2(-x, -y); }
};

v2 operator+(v2 a, v2 b) { return v2(a.x + b.x, a.y + b.y); }
v2 operator-(v2 a, v2 b) { return v2(a.x - b.x, a.y - b.y); }
v2 operator*(v2 a, float b) { return v2(a.x * b, a.y * b); }
float dot(v2 a, v2 b) { return a.x * b.x + a.y * b.y; }
v2 norm(v2 n) { float l = sqrtf(dot(n, n)); return n * (1.0f / l); }
v2 c2CCW90(v2 a) { return v2(a.y, -a.x); }

struct halfspace
{
	halfspace() { }
	halfspace(v2 n, float d) { this->n = n; this->d = d; }
	v2 n;
	float d;
};

float distance(halfspace h, v2 p) { return dot(h.n, p) - h.d; }
v2 intersect(v2 a, v2 b, float da, float db) { return a + (b - a) * (da / (da - db)); }
v2 intersect(halfspace h, v2 a, v2 b) { return intersect(a, b, distance(h, a), distance(h, b)); }

using poly = std::vector<v2>;

// -------------------------------------------------------------------------------------------------
// PNG saving functions from Richard Mitton's tigr (tiny graphics) library.

#include <stdio.h>

struct cp_pixel_t
{
	uint8_t r;
	uint8_t g;
	uint8_t b;
	uint8_t a;
};

struct cp_image_t
{
	int w;
	int h;
	cp_pixel_t* pix;
};

static uint8_t cp_paeth(uint8_t a, uint8_t b, uint8_t c)
{
	int p = a + b - c;
	int pa = abs(p - a);
	int pb = abs(p - b);
	int pc = abs(p - c);
	return (pa <= pb && pa <= pc) ? a : (pb <= pc) ? b : c;
}

typedef struct cp_save_png_data_t
{
	uint32_t crc;
	uint32_t adler;
	uint32_t bits;
	uint32_t prev;
	uint32_t runlen;
	FILE *fp;
} cp_save_png_data_t;

uint32_t CP_CRC_TABLE[] = {
	0, 0x1db71064, 0x3b6e20c8, 0x26d930ac, 0x76dc4190, 0x6b6b51f4, 0x4db26158, 0x5005713c,
	0xedb88320, 0xf00f9344, 0xd6d6a3e8, 0xcb61b38c, 0x9b64c2b0, 0x86d3d2d4, 0xa00ae278, 0xbdbdf21c
};

static void cp_put8(cp_save_png_data_t* s, uint32_t a)
{
	fputc(a, s->fp);
	s->crc = (s->crc >> 4) ^ CP_CRC_TABLE[(s->crc & 15) ^ (a & 15)];
	s->crc = (s->crc >> 4) ^ CP_CRC_TABLE[(s->crc & 15) ^ (a >> 4)];
}

static void cp_update_adler(cp_save_png_data_t* s, uint32_t v)
{
	uint32_t s1 = s->adler & 0xFFFF;
	uint32_t s2 = (s->adler >> 16) & 0xFFFF;
	s1 = (s1 + v) % 65521;
	s2 = (s2 + s1) % 65521;
	s->adler = (s2 << 16) + s1;
}

static void cp_put32(cp_save_png_data_t* s, uint32_t v)
{
	cp_put8(s, (v >> 24) & 0xFF);
	cp_put8(s, (v >> 16) & 0xFF);
	cp_put8(s, (v >> 8) & 0xFF);
	cp_put8(s, v & 0xFF);
}

static void cp_put_bits(cp_save_png_data_t* s, uint32_t data, uint32_t bitcount)
{
	while (bitcount--)
	{
		uint32_t prev = s->bits;
		s->bits = (s->bits >> 1) | ((data & 1) << 7);
		data >>= 1;

		if (prev & 1)
		{
			cp_put8(s, s->bits);
			s->bits = 0x80;
		}
	}
}

static void cp_put_bitsr(cp_save_png_data_t* s, uint32_t data, uint32_t bitcount)
{
	while (bitcount--)
		cp_put_bits(s, data >> bitcount, 1);
}

static void cp_begin_chunk(cp_save_png_data_t* s, const char* id, uint32_t len)
{
	cp_put32(s, len);
	s->crc = 0xFFFFFFFF;
	cp_put8(s, id[0]);
	cp_put8(s, id[1]);
	cp_put8(s, id[2]);
	cp_put8(s, id[3]);
}

static void cp_encode_literal(cp_save_png_data_t* s, uint32_t v)
{
	// Encode a literal/length using the built-in tables.
	// Could do better with a custom table but whatever.
	     if (v < 144) cp_put_bitsr(s, 0x030 + v -   0, 8);
	else if (v < 256) cp_put_bitsr(s, 0x190 + v - 144, 9);
	else if (v < 280) cp_put_bitsr(s, 0x000 + v - 256, 7);
	else              cp_put_bitsr(s, 0x0c0 + v - 280, 8);
}

static void cp_encode_len(cp_save_png_data_t* s, uint32_t code, uint32_t bits, uint32_t len)
{
	cp_encode_literal(s, code + (len >> bits));
	cp_put_bits(s, len, bits);
	cp_put_bits(s, 0, 5);
}

static void cp_end_run(cp_save_png_data_t* s)
{
	s->runlen--;
	cp_encode_literal(s, s->prev);

	if (s->runlen >= 67) cp_encode_len(s, 277, 4, s->runlen - 67);
	else if (s->runlen >= 35) cp_encode_len(s, 273, 3, s->runlen - 35);
	else if (s->runlen >= 19) cp_encode_len(s, 269, 2, s->runlen - 19);
	else if (s->runlen >= 11) cp_encode_len(s, 265, 1, s->runlen - 11);
	else if (s->runlen >= 3) cp_encode_len(s, 257, 0, s->runlen - 3);
	else while (s->runlen--) cp_encode_literal(s, s->prev);
}

static void cp_encode_byte(cp_save_png_data_t *s, uint8_t v)
{
	cp_update_adler(s, v);

	// Simple RLE compression. We could do better by doing a search
	// to find matches, but this works pretty well TBH.
	if (s->prev == v && s->runlen < 115) s->runlen++;

	else
	{
		if (s->runlen) cp_end_run(s);

		s->prev = v;
		s->runlen = 1;
	}
}

static void cp_save_header(cp_save_png_data_t* s, cp_image_t* img)
{
	fwrite("\211PNG\r\n\032\n", 8, 1, s->fp);
	cp_begin_chunk(s, "IHDR", 13);
	cp_put32(s, img->w);
	cp_put32(s, img->h);
	cp_put8(s, 8); // bit depth
	cp_put8(s, 6); // RGBA
	cp_put8(s, 0); // compression (deflate)
	cp_put8(s, 0); // filter (standard)
	cp_put8(s, 0); // interlace off
	cp_put32(s, ~s->crc);
}

static cp_pixel_t cp_make_pixel_a(uint8_t r, uint8_t g, uint8_t b, uint8_t a)
{
	cp_pixel_t p;
	p.r = r; p.g = g; p.b = b; p.a = a;
	return p;
}

static cp_pixel_t cp_make_pixel(uint8_t r, uint8_t g, uint8_t b)
{
	cp_pixel_t p;
	p.r = r; p.g = g; p.b = b; p.a = 0xFF;
	return p;
}

static long cp_save_data(cp_save_png_data_t* s, cp_image_t* img, long dataPos)
{
	cp_begin_chunk(s, "IDAT", 0);
	cp_put8(s, 0x08); // zlib compression method
	cp_put8(s, 0x1D); // zlib compression flags
	cp_put_bits(s, 3, 3); // zlib last block + fixed dictionary

	for (int y = 0; y < img->h; ++y)
	{
		cp_pixel_t *row = &img->pix[y * img->w];
		cp_pixel_t prev = cp_make_pixel_a(0, 0, 0, 0);

		cp_encode_byte(s, 1); // sub filter
		for (int x = 0; x < img->w; ++x)
		{
			cp_encode_byte(s, row[x].r - prev.r);
			cp_encode_byte(s, row[x].g - prev.g);
			cp_encode_byte(s, row[x].b - prev.b);
			cp_encode_byte(s, row[x].a - prev.a);
			prev = row[x];
		}
	}

	cp_end_run(s);
	cp_encode_literal(s, 256); // terminator
	while (s->bits != 0x80) cp_put_bits(s, 0, 1);
	cp_put32(s, s->adler);
	long dataSize = (ftell(s->fp) - dataPos) - 8;
	cp_put32(s, ~s->crc);

	return dataSize;
}

int cp_save_png(const char* file_name, const cp_image_t* img)
{
	cp_save_png_data_t s;
	long dataPos, dataSize, err;

	FILE* fp = fopen(file_name, "wb");
	if (!fp) return 1;

	s.fp = fp;
	s.adler = 1;
	s.bits = 0x80;
	s.prev = 0xFFFF;
	s.runlen = 0;

	cp_save_header(&s, (cp_image_t*)img);
	dataPos = ftell(s.fp);
	dataSize = cp_save_data(&s, (cp_image_t*)img, dataPos);

	// End chunk.
	cp_begin_chunk(&s, "IEND", 0);
	cp_put32(&s, ~s.crc);

	// Write back payload size.
	fseek(fp, dataPos, SEEK_SET);
	cp_put32(&s, dataSize);

	err = ferror(fp);
	fclose(fp);
	return !err;
}

// -------------------------------------------------------------------------------------------------
// Raster functions Richard Mitton's tigr (tiny graphics) library.

#include <math.h>

using pixel = cp_pixel_t;
using img = cp_image_t;

pixel color_red() { return { 0xFF, 0, 0, 0xFF }; }
pixel color_green() { return { 0, 0xFF, 0, 0xFF }; }
pixel color_blue() { return { 0, 0, 0xFF, 0xFF }; }
pixel color_white() { return { 0xFF, 0xFF, 0xFF, 0xFF }; }
pixel color_black() { return { 0, 0, 0, 0xFF }; }
pixel color_clear() { return { 0, 0, 0, 0 }; }

img blank(int w, int h)
{
	img bmp;
	bmp.w = w;
	bmp.h = h;
	bmp.pix = (pixel*)malloc(sizeof(pixel) * w * h);
	return bmp;
}

int save(const char* file_name, const img* bmp)
{
	return cp_save_png(file_name, (const img*)bmp);
}

// Expands 0-255 into 0-256
#define EXPAND(X) ((X) + ((X) > 0))

void plot(img* bmp, int x, int y, pixel pix)
{
	if (x >= 0 && y >= 0 && x < bmp->w && y < bmp->h) {
		int xa = EXPAND(pix.a);
		int a = xa * xa;
		int i = y * bmp->w + x;
		bmp->pix[i].r += (unsigned char)((pix.r - bmp->pix[i].r) * a >> 16);
		bmp->pix[i].g += (unsigned char)((pix.g - bmp->pix[i].g) * a >> 16);
		bmp->pix[i].b += (unsigned char)((pix.b - bmp->pix[i].b) * a >> 16);
		bmp->pix[i].a += (unsigned char)((pix.a - bmp->pix[i].a) * a >> 16);
	}
}

void clear(img* bmp, pixel color)
{
	int count = bmp->w * bmp->h;
	int n;
	for (n = 0; n < count; n++) {
		bmp->pix[n] = color;
	}
}

void line(img* bmp, int x0, int y0, int x1, int y1, pixel color)
{
	int sx, sy, dx, dy, err, e2;
	dx = abs(x1 - x0);
	dy = abs(y1 - y0);
	if (x0 < x1) sx = 1;
	else sx = -1;
	if (y0 < y1) sy = 1;
	else sy = -1;
	err = dx - dy;

	do {
		plot(bmp, x0, y0, color);
		e2 = 2 * err;
		if (e2 > -dy) {
			err -= dy;
			x0 += sx;
		}
		if (e2 < dx) {
			err += dx;
			y0 += sy;
		}
	} while ((x0 != x1) | (y0 != y1));
}

void circle(img* bmp, int x0, int y0, int r, pixel color) {
	int E = 1 - r;
	int dx = 0;
	int dy = -2 * r;
	int x = 0;
	int y = r;

	plot(bmp, x0, y0 + r, color);
	plot(bmp, x0, y0 - r, color);
	plot(bmp, x0 + r, y0, color);
	plot(bmp, x0 - r, y0, color);

	while (x < y - 1) {
		x++;

		if (E >= 0) {
			y--;
			dy += 2;
			E += dy;
		}

		dx += 2;
		E += dx + 1;

		plot(bmp, x0 + x, y0 + y, color);
		plot(bmp, x0 - x, y0 + y, color);
		plot(bmp, x0 + x, y0 - y, color);
		plot(bmp, x0 - x, y0 - y, color);

		if (x != y) {
			plot(bmp, x0 + y, y0 + x, color);
			plot(bmp, x0 - y, y0 + x, color);
			plot(bmp, x0 + y, y0 - x, color);
			plot(bmp, x0 - y, y0 - x, color);
		}
	}
}

// -------------------------------------------------------------------------------------------------
// Demo program to split a shape and render a png image of the output.

#include <assert.h>

// Width/height of the output image.
int w = 128;
int h = 128;
img* vis; // Output image.

// World to output image x.
int w2x(v2 p)
{
	int x = (int)p.x + w / 2;
	return x;
}

// World to output image y.
int w2y(v2 p)
{
	int y = (int)p.y + h / 2;
	return h - y;
}

void draw_poly(img* bmp, poly p, pixel color)
{
	v2 a = p[p.size() - 1];
	int ax = w2x(a);
	int ay = w2y(a);
	for (int i = 0; i < (int)p.size(); ++i) {
		v2 b = p[i];
		int bx = w2x(b);
		int by = w2y(b);
		line(bmp, ax, ay, bx, by, color);
		a = b;
		ax = bx;
		ay = by;
	}
}

struct sutherland_hodgman_output
{
	poly front;
	poly back;
	std::vector<v2> intersections;
};

bool in_front(float distance, float epsilon) { return distance > epsilon; }
bool behind(float distance, float epsilon) { return distance < -epsilon; }
bool on(float distance, float epsilon) { return !in_front(distance, epsilon) && !behind(distance, epsilon); }

// See: https://gamedevelopment.tutsplus.com/tutorials/how-to-dynamically-slice-a-convex-shape--gamedev-14479
sutherland_hodgman_output sutherland_hodgman(halfspace split, poly in, const float k_epsilon = 1.e-4f)
{
	sutherland_hodgman_output out;
	v2 a = in.back();
	float da = distance(split, a);

	for(int i = 0; i < (int)in.size(); ++i) {
		v2 b = in[i];
		float db = distance(split, b);

		if(in_front(db, k_epsilon)) {
			if(behind(da, k_epsilon)) {
				v2 i = intersect(b, a, db, da);
				out.front.push_back(i);
				out.back.push_back(i);
				out.intersections.push_back(i);
			}
			out.front.push_back(b);
		} else if(behind(db, k_epsilon)) {
			if(in_front(da, k_epsilon)) {
				v2 i = intersect(a, b, da, db);
				out.front.push_back(i);
				out.back.push_back(i);
				out.intersections.push_back(i);
			} else if(on(da, k_epsilon)) {
				out.back.push_back(a);
			}
			out.back.push_back(b);
		} else {
			out.front.push_back(b);
			if(on(da, k_epsilon)) {
				out.back.push_back(b);
			}
		}

		a = b;
		da = db;
	}

	return out;
}

std::vector<v2> clip(poly p0, poly p1)
{
	v2 a = p0[p0.size() - 1];
	for (int i = 0; i < (int)p0.size(); ++i) {
		v2 b = p0[i];
		v2 n = -norm(c2CCW90(b - a));
		halfspace h = halfspace(n, dot(n, a));
		sutherland_hodgman_output sho = sutherland_hodgman(h, p1);
		p1 = sho.back;
		a = b;
	}
	return p1;
}

bool near_surface(poly a, v2 b, const float k_epsilon = 1.0e-4f)
{
	v2 v0 = a[a.size() - 1];
	for (int i = 0; i < (int)a.size(); ++i) {
		v2 v1 = a[i];
		v2 n = norm(c2CCW90(v1 - v0));
		halfspace h = halfspace(n, dot(n, v0));
		float d = distance(h, b);
		if (abs(d) < k_epsilon) return true;
		v0 = v1;
	}
	return false;
}

std::vector<v2> intersection_points(poly p0, poly p1)
{
	std::vector<v2> clipped = clip(p0, p1);
	std::vector<v2> out;
	for (int i = 0; i < (int)clipped.size(); ++i) {
		v2 p = clipped[i];
		if (near_surface(p0, p) && near_surface(p1, p)) {
			out.push_back(p);
		}
	}
	return out;
}

int main()
{
	poly a = {
		{ -20.0f, -20.0f },
		{ -60.0f, 0.0f },
		{ 40.0f, 20.0f },
		{ 50.0f, -10.0f },
	};

	poly b = {
		{ -30.0f, 45.0f },
		{ 15, 30.0f },
		{ -50, -30.0f },
	};

	img bmp = blank(w, h);
	vis = &bmp;

	clear(vis, color_white());
	draw_poly(vis, a, color_blue());
	draw_poly(vis, b, color_red());

	std::vector<v2> out = intersection_points(a, b);
	for (int i = 0; i < (int)out.size(); ++i) {
		circle(vis, w2x(out[i]), w2y(out[i]), 5, color_green());
	}

	save("out.png", vis);

	free(vis->pix);

	return 0;
}

That's pretty cool! Thanks for sharing. By the way, shouldn't the filename be something like edge_intersections.cpp since this is C++?

Yes it’s c++

Here are some renders