marmakoide’s gists

marmakoide / complex-step-jacobian.py

Created August 19, 2020 12:55

Computes the jacobian of a real-valued function, using the complex step method. The result is accurate nearly down to the machine epsilon with a very simple code, at the expense of uncessary computations, and assuming that the function to differentiate is defined in the complex domain

	def get_jacobian(func, X, h = 1e-7):
	J = []
	for i in range(X.shape[0]):
	Xp = X.astype(complex)
	Xp[i] += h * 1j
	J.append([numpy.imag(func(Xp))])

	J = numpy.concatenate(J, axis = 0)
	J = numpy.transpose(J)
	J /= h

marmakoide / bevel-box.scad

Last active June 12, 2020 08:55

Bevel rectangle and boxes for OpenSCAD

	module bevel_rectangle(size, bevel) {
	s = (size[0] + size[1] - 2 * bevel) / sqrt(2);

	intersection() {
	square(size, center = true);
	rotate(45)
	square([s, s], center = true);
	}
	}

marmakoide / reverse-euclidean-distance-transform.py

Created May 21, 2020 09:13

SIMD implementation of 2d reverse euclidean distance. See "Optimal Separable Algorithms to Compute the Reverse Euclidean Distance Transformation and Discrete Medial Axis in Arbitrary Dimension" D. Coeurjolly et al.

	import numpy



	def reverse_euclidean_distance_transform(F):
	ret = numpy.empty_like(F, dtype = int)

	def get_parabola_matrix(N):
	X = numpy.arange(N, dtype = int)
	G = numpy.empty((N, N), dtype = int)

marmakoide / euclidean-distance-transform-2d.py

Created May 21, 2020 09:09

SIMD implementation of 2d euclidean distance. For each non-zero input pixel, returns the exact squared distance to the closest zero input pixel.

	import numpy

	def euclidean_distance_transform(F):
	ret = numpy.empty_like(F, dtype = int)

	def get_parabola_matrix(N):
	X = numpy.arange(N, dtype = int)
	G = numpy.empty((N, N), dtype = int)
	G[:] = X
	G -= X[:, None]

marmakoide / oric-lores-color-attributes.bas

Created April 11, 2020 20:24

Demonstration of color attributes usage in low resolution mode, Oric 1 computer

	10 FOR I=1 TO 27
	20 POKE #BB80+I*40,16
	30 POKE #BB81+I*40,I AND 7
	40 FOR J=0 TO 38
	50 POKE #BB82+I*40+J,126
	60 NEXT J
	70 POKE #BB80+I*40+20,23
	80 NEXT I

marmakoide / dft-2d-rot.py

Last active March 31, 2020 09:19

Rotation of 2d DFT coefficients, avoiding computing a DFT by directly transforming the DFT coefficients

	import numpy

	def rot180_fft_2d(X, N):
	I = numpy.repeat([numpy.arange(N)], N, axis = 0)
	K = numpy.exp((2 * numpy.pi / N) * (I + I.T) * 1j)
	return K * numpy.conj(X)

	def rot180_rfft_2d(X, N):
	I = numpy.repeat([numpy.arange(N)], N, axis = 0)
	K = numpy.exp((2 * numpy.pi / N) * (I + I.T) * 1j)[:,:N // 2 + 1]

marmakoide / axis-angle-rotation-matrix.py

Last active December 13, 2022 13:21

Computes the rotation matrix equivalent to an axis & angle rotation

	import numpy

	def axis_angle_to_rotation_matrix(U, theta):
	cos_theta, sin_theta = numpy.cos(theta, dtype = U.dtype), numpy.sin(theta, dtype = U.dtype)
	U_x = numpy.array([[0, -U[2], U[1]], [U[2], 0, -U[0]], [-U[1], U[0], 0]], dtype = U.dtype)
	return cos_theta * numpy.identity(3, dtype = U.dtype) + sin_theta * U_x + (1 - cos_theta) * numpy.outer(U, U)

marmakoide / hexacosichoron.py

Created February 12, 2020 07:27

Generates the vertices of a unit hexacosichoron

	import itertools
	from sympy.combinatorics import Permutation

	def get_hexacosichoron_vertices():
	phi = (1 + 5 ** .5) / 2
	U = [.5 * phi, .5, .5 / phi]

	permutation_list = [P for P in itertools.permutations(range(4)) if Permutation(P).is_even]

	for i in range(16):

marmakoide / bridson_algorithm.py

Created November 19, 2019 20:37

An implementation of Bridson's algorithm for blue noise

	import numpy

	def sample_annulus(r1, r2, k):
	theta = numpy.random.rand(k)
	theta = 2 numpy.pi
	R = numpy.sqrt(numpy.random.rand(k))
	R *= r2 - r1
	R += r1
	S = numpy.array([numpy.cos(theta), numpy.sin(theta)]).T
	S *= R[:,None]

marmakoide / quickhull-2d.py

Created January 17, 2019 07:12

A nice, short compact 2d convex hull routine, implementing the Quickhull algorithm

	import numpy



	def quickhull_2d(S):
	def process(P, a, b):
	signed_dist = numpy.cross(S[P] - S[a], S[b] - S[a])
	K = [i for s, i in zip(signed_dist, P) if s > 0 and i != a and i != b]

	if len(K) == 0:

Devert Alexandre marmakoide