spiral-3d.py

### Copyright (c) 2015, Sergey Shishmintzev
### License: Apache/GPLv3
### Example video:  http://youtu.be/aJ9IRxgpH10
### Interpreter: Python 2.7.5 with NumPy 1.9.1 and Pillow 2.7.0

power = 19

name = "spiral"

frame_count = 15000

scale = 2

image_width = 1080
image_height = 1080

alpha_z = 55

import numpy as np
import numpy.random as nr
from PIL import Image
import colorsys

np.set_printoptions(linewidth = 180, edgeitems=10, suppress = True)

q_zero = (0, 0, 0, 0)
q_unit = (1, 0, 0, 0)
q_i = (0, 1, 0, 0)
q_j = (0, 0, 1, 0)
q_k = (0, 0, 0, 1)

def q_mult(q1, q2):
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    w = w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2
    x = w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2
    y = w1 * y2 + y1 * w2 + z1 * x2 - x1 * z2
    z = w1 * z2 + z1 * w2 + x1 * y2 - y1 * x2
    return w, x, y, z

def q_abs(q):
  return np.sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3])

def q_exp(q):
    w, x, y, z = q
    modulus = np.exp(w)
    vector_modulus = np.sqrt(x*x + y*y + z*z)
    #K = modulus*np.sin(vector_modulus)/vector_modulus
    K = modulus*np.sinc(vector_modulus/np.pi)
    return modulus*np.cos(vector_modulus), K*x, K*y, K*z

def additive_binary_construction(components, n):
    result = 0
    bit_index = 0
    while n > 0:
        if (n & 1) == 1:
            result += components[bit_index]
        n >>= 1
        bit_index += 1
    return result

def multiplicative_binary_construction(components, n):
    result = 1
    bit_index = 0
    while n > 0:
        if (n & 1) == 1:
              result *= components[bit_index]
        n >>= 1
        bit_index += 1
    return result

# generate palette
palette = [0] * 256
for i, h in enumerate(np.linspace(0, 3, 256)):
    c = colorsys.hsv_to_rgb(h, 1.0, 1.0)
    palette[i] = (int(c[0]*255), int(c[1]*255), int(c[2]*255))

def alphablend(dst, src, alpha):
    return ((dst[0]*(255-alpha) + src[0]*alpha)>>8, (dst[1]*(255-alpha) + src[1]*alpha)>>8, (dst[2]*(255-alpha) + src[2]*alpha)>>8)


r = np.linspace(1.0/power, 1.0, power)
# np.exp(0.5*np.log(1+r[i]**2) + 1j*np.log(r[i]+1j).imag) == r[i] + 1j
spiral_components = [q_exp((0.5*np.log(1+r[i]**2), np.log(r[i]+1j).imag, 0.34657359028, 0)) for i in xrange(power)]


print "spiral_components:"
print spiral_components
modulus = [q_abs(q) for q in spiral_components]
print "component modulus:"
print modulus

spiral_modulus = np.prod(modulus)
print "spiral_modulus =", spiral_modulus

hamming_weight_components = np.ones(power)

# Recombination examples
#recombination_components = nr.randint(0, 2, power)
#recombination_components = np.ones(power)
#recombination_components = np.arange(power)+1
#recombination_components = np.array([1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3,])
#recombination_components = np.array([1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1,])
#recombination_components = np.array([1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5,])
recombination_components = np.array([1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10,])
#recombination_components = np.array([5, 4, 4, 3, 3, 2, 2, 1, 1, 0, -1, -1, -2, -2, -3, -3, -4, -4, -5])

assert len(recombination_components) == power, len(recombination_components)
print "recombination_components", recombination_components

print "Generate data..."


spiral_points = np.empty((2**power, 4), np.float64)
hamming_weight = np.empty(2**power, np.uint16)
recombination = np.empty(2**power, np.int32)

# slow code

"""
for i in xrange(2**power):
    spiral_points[i] = multiplicative_binary_construction(spiral_components, i)
    hamming_weight[i] = additive_binary_construction(hamming_weight_components, i)
    recombination[i] = additive_binary_construction(recombination_components, i)
"""

# fast code (based boolean cube walking)

spiral_points[0] = q_unit
hamming_weight[0] = 0
recombination[0] = 0

pair_step = 2**(power-1)
for i in xrange(power):
    left_pos = 0
    spiral_component = spiral_components[-i-1]
    hamming_weight_component = hamming_weight_components[-i-1]
    recombination_component = recombination_components[-i-1]
    for j in xrange(2**i):
        spiral_points[left_pos+pair_step] = q_mult(spiral_points[left_pos], spiral_component)
        hamming_weight[left_pos+pair_step] = hamming_weight[left_pos] + hamming_weight_component
        recombination[left_pos+pair_step] = recombination[left_pos] + recombination_component
        left_pos += 2*pair_step
    pair_step /= 2


rotation = q_exp((0.0, 1.0/(frame_count/10.0), 1.0/(frame_count/10.0), 4*np.pi/(frame_count/10.0)))

recombination_widdle_factor = np.empty((2**power, 4), dtype=np.float64)
for i in xrange(2**power):
  recombination_widdle_factor[i,:] = q_mult(rotation, q_exp((0.0, 0, 2*np.pi*recombination[i]/frame_count, 0)))

print "Render..."

step_x = 1.0*(image_width-1) / (2*spiral_modulus)
step_y = 1.0*(image_height-1) / (2*spiral_modulus)

color_index = np.uint16(255*hamming_weight/power)

qcompoment = 1

for frame in xrange(frame_count+1):
    print "frame #", frame
    image = Image.new("RGB", (image_width, image_height), (255, 255, 255))
    pixels = image.load()

    spiral_points_x = np.int32((scale*spiral_points[:,0]+spiral_modulus)*step_x)
    spiral_points_x = np.clip(spiral_points_x, 0, image_width-1)
    spiral_points_y = np.int32((spiral_modulus-scale*spiral_points[:,qcompoment])*step_y)
    spiral_points_y = np.clip(spiral_points_y, 0, image_height-1)

    # render
    for i in xrange(2**power):
        x = spiral_points_x[i]
        y = spiral_points_y[i]
        pixels[x, y] = alphablend(pixels[x, y], palette[color_index[i]], alpha_z)

    image.save("%s-frame%05d-z.png"%(name, frame), "PNG")

    for i in xrange(2**power):
      spiral_points[i] = q_mult(spiral_points[i], recombination_widdle_factor[i])