# Implementation of Spiral Architecture¶

Spiral architecture has been introduced by Sheridan in

The implementation with hyperpels that we use in the following is presented in

## Inizializations and auxiliary functions¶

In [2]:
import numpy as np
from numpy import fft
from numpy import linalg as LA
from scipy import ndimage
# from scipy import linalg
from scipy import signal
import matplotlib.pyplot as plt
import matplotlib.cm as cm

import os

In [3]:
%matplotlib inline

In [4]:
# ALPHABET = "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"
ALPHABET = "0123456"

def base7_encode(num, alphabet=ALPHABET):
"""Encode a number in Base X

num: The number to encode
alphabet: The alphabet to use for encoding
"""
if (num == 0):
return alphabet[0]
arr = []
base = len(alphabet)
while num:
rem = num % base
num = num // base
arr.append(alphabet[rem])
arr.reverse()
return int(''.join(arr))

def base7_decode(string, alphabet=ALPHABET):
"""Decode a Base X encoded string into the number

Arguments:
- string: The encoded string
- alphabet: The alphabet to use for encoding
"""
base = len(alphabet)
strlen = len(string)
num = 0

idx = 0
for char in string:
power = (strlen - (idx + 1))
num += alphabet.index(char) * (base ** power)
idx += 1

return num

In [5]:
def rgb2gray(rgb):
r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
return gray

In [6]:
example =  1 - rgb2gray(plt.imread('./test-images/butterfly.png'))
plt.imshow(example,cmap=cm.binary)

Out[6]:
<matplotlib.image.AxesImage at 0x10a86d350>

## Oversampling¶

We start with a $M\times M$ image, and transform it to a $7M \times 7M$ one

In [7]:
def oversampling(image, factor = 7):
old_shape = image.shape
new_shape = (factor*old_shape[0], factor*old_shape[1])
new_image = np.zeros(new_shape)
for i in range(old_shape[0]):
for j in range(old_shape[1]):
new_image[factor*i:factor*i+factor,factor*j:factor*j+factor] = image[i,j]*np.ones((factor,factor))
return new_image

In [8]:
oversampled_example = oversampling(example)
plt.imshow(oversampled_example,cmap = cm.binary)

Out[8]:
<matplotlib.image.AxesImage at 0x10cbec350>

## Rows/Columns hexagonal grid¶

The grid where we place the hyperpel is centered in the center of the central hyperpel.

In [9]:
# Distance between centers of adjacent hexagons
# with the hyperpel
#
# o o x x x x x o o
# o x x x x x x x o
# o x x x x x x x o
# x x x x x x x x x
# x x x x x x x x x
# o x x x x x x x o
# o x x x x x x x o
# o o x x x x x o o

dist = np.sqrt(65)

# The centered hyperpel
hyperpel = np.array([\
[-1,4],[0,4],[1,4],[2,4],[3,4],\
[-2,3],[-1,3], [0,3], [1,3], [2,3], [3,3], [4,3],\
[-2,2],[-1,2], [0,2], [1,2], [2,2], [3,2], [4,2],\
[-3,1],[-2,1],[-1,1], [0,1], [1,1], [2,1], [3,1], [4,1],[5,1],\
[-3,0],[-2,0],[-1,0], [0,0], [1,0], [2,0], [3,0], [4,0],[5,0],\
[-2,-1],[-1,-1], [0,-1], [1,-1], [2,-1], [3,-1], [4,-1],\
[-2,-2],[-1,-2], [0,-2], [1,-2], [2,-2], [3,-2], [4,-2],\
[-1,-3], [0,-3], [1,-3], [2,-3], [3,-3]])

for i in range(56):
hyperpel[i] = [hyperpel[i,1],hyperpel[i,0]]

In [10]:
# This function returns the coordinates of an hexagon, given
# its position with row/column, for the above hyperpel
# and assuming an oversampled image by a factor of 7

# Here the hex of position [0,1] is right UP w.r.t. the central hexagon

def hex2cartesian(row, column):
x = 7 * column
# even column
if column % 2 == 0:
y = 8*row
# odd row
else:
y = 4 + 8*row

return x,y

In [11]:
# This function returns the coordinates of the points of
# the hyperpel corresponding to the hexagon of given
# row/column positions

# Here, we assume the center of the [0,0] hyperpel to be at the
# intersection of the 2 rightmost pixels of the central 4,
# while the "center of the image" will be as in picture:
# o o x x x x x o o
# o x x x x x x x o
# o x x x x x x x o
# x x x x x x x x x
# x x x C x x x x x
# o x x x x x x x o
# o x x x x x x x o
# o o x x x x x o o

def center_hyperpel(image,column,row):
# The center of the image are 4 pixels, we denote as center
# the one left-down:
#
#  O O
#  C O
cent_y = int(np.floor(float(image.shape[0])/2))
cent_x = int(np.floor(float(image.shape[1])/2))

h_x,h_y = hex2cartesian(row,column)

# We return the translated one
return hyperpel + [cent_y,cent_x] + [h_y,h_x]


The following shows that indeed the hyperpels are generated correctly

In [12]:
ciao = np.ones((40,40))
hp = center_hyperpel(ciao,0,0)
for i in range(56):
ciao[hp[i,0],hp[i,1]] =0
hp = center_hyperpel(ciao,1,0)
for i in range(56):
ciao[hp[i,0],hp[i,1]] =0
hp = center_hyperpel(ciao,2,0)
for i in range(56):
ciao[hp[i,0],hp[i,1]] =0
plt.imshow(ciao,cmap=cm.binary)

Out[12]:
<matplotlib.image.AxesImage at 0x11f9514d0>
In [13]:
# This function returns the value of the hexagonal pixel,
# depending on the values of the oversampled image
# given its row/column position, for the above hyperpel

def hex_value(oversampled_image,row,column):
hp = center_hyperpel(oversampled_image,row,column)
val = 0.
for i in range(56):
val += oversampled_image[hp[i,0],hp[i,1]]

return val/56

In [14]:
# This function returns the grid for the hex row/column representation
# given the ORIGINAL size of the image

def hex_grid(image_size = 256):
up_rows = down_rows = int(np.floor((7*image_size/2-4)/8))
left_columns = int(np.floor((7*image_size/2-4)/9))
right_columns = int(np.floor((7*image_size/2-5)/9))

# My computations suggested +1, but up to
# +22 works without catching any exception
return np.zeros((up_rows+down_rows+25, left_columns+right_columns+25))

In [15]:
# This function, given an image, fills the corresponding hex grid

def image2hex(image):
original_size = image.shape[0]
oversampled_image = oversampling(image)
grid = hex_grid(original_size)

hex_center = [(grid.shape[0]-1)/2,(grid.shape[1]-1)/2]

for i in range(grid.shape[0]):
for j in range(grid.shape[1]):
try:
grid[i,j] = hex_value(oversampled_image,i-hex_center[0], j-hex_center[1])
except IndexError:
continue

return grid

In [16]:
def hex2image(grid, image_size=256):
image = np.zeros((image_size*7,image_size*7))
hex_center = [(grid.shape[0]-1)/2,(grid.shape[1]-1)/2]
for i in range(grid.shape[0]):
for j in range(grid.shape[1]):
hp = center_hyperpel(image,i-hex_center[0],j-hex_center[1])
for k in range(56):
try:
image[hp[k,0],hp[k,1]] = grid[i,j]
except IndexError:
#                     print("[{},{}]").format(hp[k,0],hp[k,1])
continue
return image

In [17]:
example_hex = hex2image(image2hex(example))
# plt.imsave("tryout.png", example_hex, cmap=cm.binary)
plt.imshow(example_hex, cmap=cm.binary)

Out[17]:
<matplotlib.image.AxesImage at 0x11fa9c0d0>

## Spiral Architecture Hexagonal grid¶

Algorithm inspired from "A New Simulation of Spiral Architecture"

In [18]:
# They use a different hyperpel!

# The centered hyperpel
hyperpel = np.array([\
[-1,4],[0,4],[1,4],[2,4],[3,4],\
[-2,3],[-1,3], [0,3], [1,3], [2,3], [3,3], [4,3],\
[-2,2],[-1,2], [0,2], [1,2], [2,2], [3,2], [4,2],\
[-3,1],[-2,1],[-1,1], [0,1], [1,1], [2,1], [3,1], [4,1],[5,1],\
[-3,0],[-2,0],[-1,0], [0,0], [1,0], [2,0], [3,0], [4,0],[5,0],\
[-2,-1],[-1,-1], [0,-1], [1,-1], [2,-1], [3,-1], [4,-1],\
[-2,-2],[-1,-2], [0,-2], [1,-2], [2,-2], [3,-2], [4,-2],\
[-1,-3], [0,-3], [1,-3], [2,-3], [3,-3]])

hyperpel_sa = hyperpel - np.array([1,1])

In [19]:
# Here, given a spiral address, we return the corresponding hexagon
# in the form row/column

# Split the number in basic unit and call the auxiliary function
# Here we reverse the order, so that the index corresponds to the
# decimal position

for i in range(len(digits)):
if int(digits[i])<0 or int(digits[i])>6:
return
elif digits[i]!= '0':

# This computes the row/column positions of the base cases,
# that is, in the form a*10^(zeros).
def sa2hex_aux(a, zeros):
# Base cases
if zeros == 0:
if a == 0:
return np.array([0,0])
elif a == 1:
return np.array([0,8])
elif a == 2:
return np.array([-7,4])
elif a == 3:
return np.array([-7,-4])
elif a == 4:
return np.array([0,-8])
elif a == 5:
return np.array([7,-4])
elif a == 6:
return np.array([7,4])

return sa2hex_aux(a,zeros-1)+ 2*sa2hex_aux(a%6 +1,zeros-1)

In [20]:
# Computes the value of the hyperpel corresponding to the given
# spiral coordinate.

# WORKS ONLY FOR IMAGES WITH EVEN PIXELS
center = [(oversampled_image.shape[0])/2,(oversampled_image.shape[1])/2]

# The center of the hyperpel for them is the point [1,1] for us
hp = hyperpel_sa + sa2hex(spiral_address) + center
val = 0.
for i in range(56):
val += oversampled_image[hp[i,0],hp[i,1]]

return val/56

In [21]:
# Assigns the given value to the hyperpel corresponding to the given
# spiral coordinate

# WORKS ONLY FOR IMAGES WITH EVEN PIXELS
center = [(oversampled_image.shape[0])/2,(oversampled_image.shape[1])/2]

hp = hyperpel_sa + sa2hex(spiral_address) + center
for i in range(56):
oversampled_image[hp[i,0],hp[i,1]] = value

return oversampled_image



In [22]:
addition_table = [
[0,1,2,3,4,5,6],
[1,63,15,2,0,6,64],
[2,15,14,26,3,0,1],
[3,2,26,25,31,4,0],
[4,0,3,31,36,42,5],
[5,6,0,4,42,41,53],
[6,64,1,0,5,53,52]
]

In [23]:
# Auxiliary function to recover a vector with the
# digits of a given integer (in inverse order)

def int2intvec(a):
digit = a%10
vec = np.array([digit],dtype=int)
a = (a-digit)/10
while a!=0:
digit = a%10
vec = np.append(vec,int(digit))
a = (a-digit)/10
return vec

In [24]:
def spiral_add(a,b,mod=0):
dig_a = int2intvec(a)
dig_b = int2intvec(b)

if (dig_a<0).any() or (dig_a>7).any() \
or (dig_b<0).any() or (dig_b>7).any():
return

if len(dig_a) == 1 and len(dig_b)==1:

if len(dig_a) < len(dig_b):
dig_a.resize(len(dig_b))
elif len(dig_b) < len(dig_a):
dig_b.resize(len(dig_a))

res = 0

for i in range(len(dig_a)):

if i == len(dig_a)-1:
else:
res += (temp%10)*(10**i)

dig_a[i+1] = str(carry_on)

if mod!=0:
return res%mod

return res

In [25]:
%%timeit

10000 loops, best of 3: 134 µs per loop

In [26]:
%%timeit
(23 + 15)%100

10000000 loops, best of 3: 22.6 ns per loop


#### Spiral multiplication¶

In [27]:
multiplication_table = [
[0,0,0,0,0,0,0],
[0,1,2,3,4,5,6],
[0,2,3,4,5,6,1],
[0,3,4,5,6,1,2],
[0,4,5,6,1,2,3],
[0,5,6,1,2,3,4],
[0,6,1,2,3,4,5],
]

In [28]:
def spiral_mult(a,b, mod=0):
dig_a = int2intvec(a)
dig_b = int2intvec(b)

if (dig_a<0).any() or (dig_a>7).any() \
or (dig_b<0).any() or (dig_b>7).any():
return

sa_mult = int(0)

for i in range(len(dig_b)):
for j in range(len(dig_a)):
temp = multiplication_table[dig_a[j]][dig_b[i]]*(10**(i+j))

if mod!=0:
return sa_mult%mod

return sa_mult


#### Tests¶

Spiral multiplication is roughly three order slower than normal multiplication

In [29]:
%%timeit
spiral_mult(435,316,100)

100 loops, best of 3: 2.02 ms per loop

In [30]:
%%timeit
435*316 & 100

10000000 loops, best of 3: 21.5 ns per loop

In [31]:
%%timeit
spiral_mult(15,10,100)

1000 loops, best of 3: 463 µs per loop

In [32]:
%%timeit
15*10 & 100

10000000 loops, best of 3: 23.7 ns per loop

In [33]:
test = np.ones((40,40))
test = sa_put_value(0,0,ciao)
test = sa_put_value(1,0,ciao)
test = sa_put_value(2,0,ciao)
# plt.imshow(np.fft.fftshift(test),cmap=cm.binary)
plt.imshow(test,cmap=cm.binary)

Out[33]:
<matplotlib.image.AxesImage at 0x122ca8fd0>
In [34]:
test = sa_put_value(10,0,np.ones((70,70)))
for i in range(2,7):
test = sa_put_value(spiral_mult(10,i),0.5,test)
plt.imshow(test,cmap=cm.binary)

Out[34]:
<matplotlib.image.AxesImage at 0x122e11750>

Time to extract the values of an orbit:

In [35]:
%%timeit
for i in range(1,7):
sa_value(oversampled_example,spiral_mult(10,i))

1000 loops, best of 3: 1.58 ms per loop


Time to put the values of an orbit

In [36]:
%%timeit
test = np.ones((60,60))
for i in range(1,7):
sa_put_value(spiral_mult(10,i),1,test)

1000 loops, best of 3: 1.54 ms per loop


$7^5=16.807$ seems to be the biggest power of $7$ contained in a 256x256 image (recall that in total the pixels are $256^2=65.536$)

It corresponds to spiral addresses of at most 5 digits

In [46]:
test = np.ones((7*256,7*256))
for i in range(7**5):
sa_put_value(base7_encode(i),0,test)
plt.imshow(test, cmap=cm.binary)

In [43]:
test = np.zeros((7*256,7*256))
for i in range(7**5):
sa_index = base7_encode(i)
sa_put_value(sa_index,sa_value(oversampled_example,sa_index) ,test)
plt.imshow(test, cmap=cm.binary)

Out[43]:
<matplotlib.image.AxesImage at 0x12f793ad0>

### Videos¶

In [45]:
img = np.ones((256,256))
for i in range(7**3):
sa_put_value(base7_encode(i),0,img)
#     plt.imshow(img, cmap=cm.binary)
plt.imsave('_tmp%05d.png'%i, img,cmap = cm.binary)

FPS = 10

os.system("rm movie.mp4")

os.system("ffmpeg -r "+str(FPS)+" -b 1800 -i _tmp%05d.png movie.mp4")
os.system("rm _tmp*.png")

Out[45]:
0