transit.py

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

power = 19

name = "transit"

frame_count = 511

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)

def linear_interpolation(start_value, stop_value, start_offset, stop_offset, offset):
    return start_value + ((offset - start_offset) / (stop_offset - start_offset) * (stop_value - start_value))

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)

power_range = np.arange(power)
power_unit = np.linspace(0, 1.0, power)

source_components = np.exp(2j*np.pi*(2**power_range)/(2**power)) # the unit circle
target_components = 1j*power_unit-1 # the hamming spiral with vanishing sum of the elements

print "transition components:"
print source_components
print target_components

source_magnitude = np.prod(np.abs(source_components))
target_magnitude = np.prod(np.abs(target_components))
magnitude = np.max([source_magnitude, target_magnitude])

print "source magnitude =", source_magnitude
print "target magnitude =", target_magnitude

print "Generating data..."

hamming_weight_components = np.ones(power)
hamming_weight = np.empty(2**power, np.uint16)
for i in xrange(2**power):
    hamming_weight[i] = additive_binary_construction(hamming_weight_components, i)

print "Rendering..."

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

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

# initialize binary constuction
transit_points = np.empty(2**power, np.complex128)
transit_points[0] = 1.0

for frame in xrange(frame_count+1):
    print "frame #", frame

    # calculate transition for components
    transit_components_mag = linear_interpolation(np.abs(source_components), np.abs(target_components), 0.0, 1.0, frame*1.0/frame_count)
    transit_components_arg = linear_interpolation(np.angle(source_components), np.angle(target_components), 0.0, 1.0, frame*1.0/frame_count)
    transit_components = transit_components_mag*np.exp(1j*transit_components_arg)

    # fast multiplicative binary construction of the sequence transit_points
    pair_step = 2**(power-1)
    for i in xrange(power):
        left_pos = 0
        transit_component = transit_components[-i-1]
        for j in xrange(2**i):
            transit_points[left_pos+pair_step] = transit_points[left_pos] * transit_component
            left_pos += 2*pair_step
        pair_step /= 2

    # draw
    image = Image.new("RGB", (image_width, image_height), (255, 255, 255))
    pixels = image.load()

    transit_points_x = np.int32((transit_points.real+magnitude)*step_x)
    transit_points_x = np.clip(transit_points_x, 0, image_width-1)
    transit_points_y = np.int32((magnitude-transit_points.imag)*step_y)
    transit_points_y = np.clip(transit_points_y, 0, image_height-1)

    # range checks
    #assert np.all(transit_points_x >= 0)
    #assert np.all(transit_points_y >= 0)
    #assert np.all(transit_points_x < image_width)
    #assert np.all(transit_points_y < image_height)

    # render
    for i in xrange(2**power):
        x = transit_points_x[i]
        y = transit_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")

print "sum (vanishing?)", np.sum(transit_points)