marmakoide

import math

def bit_reversal_sequence(N):
    bit_count = int(math.ceil(math.log2(N)))
    for i in range(2 ** bit_count):
        bit_str = bin(i)[2:].rjust(bit_count, '0')
        bit_str = ''.join(reversed(bit_str))
        j = int(bit_str, 2)
        if j < N:
            yield j

marmakoide

import numpy

'''
Unit simplex in N dimension, as a triangular matrix,
only positive coeffs.
1st vertex is the origin, not in the returned matrix

I don't know if this algorithm have been discovered
before. I find it elegant, and I believe it can be
very accurate and stable with arbitrary precision numbers

marmakoide

import numpy

def get_orthogonal_projection(U):
    M = numpy.random.randn(U.shape[1], U.shape[1])
    M[:U.shape[0]] = U
    
    Q, R = numpy.linalg.qr(M.T)
    M[U.shape[0]:] = Q.T[U.shape[0]:]
    
    Q, R = numpy.linalg.qr(M.T)

marmakoide

def aabox_triangle_3d_intersecting(box_center, box_extents, ABC):
    # Algorithm described in "Fast 3D Triangle-Box Overlap Testing" by Tomas Akenine-Moller 
    ABC = ABC - box_center
    UVW = numpy.roll(ABC, -1, axis = 0) - ABC

    # Check along triangle edges
    for K in UVW:
        P = numpy.dot(ABC, numpy.array([0, -K[2], K[1]]))
        if max(-numpy.amax(P), numpy.amin(P)) > box_extents[1] * numpy.fabs(K[2]) + box_extents[2] * numpy.fabs(K[1]):
            return False

marmakoide

def aabox_segment_3d_intersecting(box_center, box_extents, A, B):
    M = (A + B) / 2 - box_center
    AB = B - A

    if numpy.any(numpy.fabs(M) > box_extents + .5 * numpy.fabs(AB)):
        return False

    U = AB / numpy.sqrt(numpy.sum(AB ** 2))
    abs_U = numpy.fabs(U)
    abs_W = numpy.fabs(numpy.cross(U, M))

marmakoide

import numppy
from scipy.spatial import ConvexHull

def convex_hull_contains_point(hull, P):
    for E in hull.equations:
        if numpy.dot(P, E[:-1]) > -E[-1]:
            return False

    return True    

marmakoide

import numpy
from scipy.spatial import ConvexHull


def get_convex_hull_2d_centroid(hull):
    C = numpy.zeros(2)
    for i in range(hull.vertices.shape[0]):
        Pi, Pj = hull.points[hull.vertices[i - 1]], hull.points[hull.vertices[i]]
        C += numpy.cross(Pi, Pj) * (Pi + Pj)
    return C / 6.

marmakoide

Golden dot spiral on a torus

This code sample demonstrates how to generates a set of points that covers a torus, for any number of points. The trick is to use a spiral and the golden angle.

42 points	256 points	999 points

marmakoide

import numpy
from matplotlib import pyplot as plot

N = 256
C = -0.4 + 0.6j
max_epoch_count = 1024
    
# Initialization
T = numpy.linspace(-2., 2., N)
X = numpy.tile(T, (N, 1))

marmakoide

import numpy
from scipy.spatial import ConvexHull


def is_oriented_simplex(S):
    M = numpy.concatenate([S, numpy.ones((S.shape[0], 1))], axis = 1)
    return numpy.linalg.det(M) > 0


def get_power_triangulation(S, R):