Interpolates between drawing a hilbert curve and the fourier transform of a hilbert curve

fraction_fourier_hilbert.py

import math
from numpy import array, hstack, vstack
import numpy as np
from PIL import Image


# these two methods were copied from https://github.com/dmishin/fft-image-experiments
def hilbert_indices(N):
    """Genrate 2^N x 2^N integer array, filled with values from 0 to 4^N along hilbert curve"""
    m = array([[0]], dtype=np.int)
    for i in range(N):
        d = 4**i
        m1 = vstack((hstack((m.T, m.T[::-1, ::-1] + 3*d)),
                      hstack((m+d, m+2*d))
                 )
       )
        m = m1
    return m
def hilbert_binary_diagram(N):
    """Generate 2^N x 2^N array of 0 and 1, drawing Hilbert curve"""
    m = array([[1,0],
               [0,0]], dtype = np.uint8)
    for i in range(1, N):

        a = 2**i
        m1 = np.zeros((2*a, 2*a), dtype = m.dtype)

        m1[0:a, 0:a] = m.T
        m1[a:2*a, 0:a] = m
        m1[a:2*a, a:2*a] = m
        m1[0:a, a:(2*a-1) ] = m.T[:, a-2::-1]

        #put additional ones
        m1[a-1,0] = 1
        m1[a, a-1] = 1
        m1[a-1, 2*a-2] = 1

        m = m1
    return m


# this method was derived from https://algassert.com/post/1710
def fast_fractional_fourier_transform_2d(vec, exponent):
    # Compute Fourier transform powers of vec.
    f0 = np.array(vec)
    f1 = np.fft.ifft2(f0, norm='ortho')
    # The following two lines could be made O(n) by re-ordering entries, but meh
    f2 = np.fft.ifft2(f1, norm='ortho')
    f3 = np.fft.ifft2(f2, norm='ortho')

    # Derive eigenbasis vectors from vec's Fourier transform powers.
    b0 = f0 + f1 + f2 + f3
    b1 = f0 + 1j*f1 - f2 - 1j*f3
    b2 = f0 - f1 + f2 - f3
    b3 = f0 - 1j*f1 - f2 + 1j*f3
    # Note: vec == (b0 + b1 + b2 + b3) / 4

    # Phase eigenbasis vectors by their eigenvalues raised to the exponent.
    b1 *= 1j**exponent
    b2 *= 1j**(2*exponent)
    b3 *= 1j**(3*exponent)

    # Recombine transformed basis vectors to get transformed vec.
    return (b0 + b1 + b2 + b3) / 4


def draw(power_of_two_zoom: int):
    bmp = hilbert_binary_diagram(power_of_two_zoom)
    spread = np.linspace(1, 0, math.prod(bmp.shape)).reshape(bmp.shape)
    f = fast_fractional_fourier_transform_2d(bmp, spread)

    # Not really a principled complex-to-color mapping... just *something*
    r = np.real(f)
    g = -np.imag(f) - np.real(f)
    b = np.imag(f)

    def normalize(v):
        v = np.minimum(v, 1)
        v = np.maximum(v, 0)
        v = (v * 255).astype(np.uint8)
        v = v.reshape((*bmp.shape, 1))
        return v

    def combo(a, b, c):
        return np.concatenate((normalize(a), normalize(b), normalize(c)), axis=2)
    rgb = combo(r, g, b)

    # img = Image.fromarray(hsv, mode="HSV")
    # img = img.convert(mode="RGB")
    img = Image.fromarray(rgb, mode="RGB")

    img.save('tmp.bmp')
    print("wrote tmp.bmp")


if __name__ == '__main__':
    draw(power_of_two_zoom=12)