Escher Red Ants¶
‘Making hard things possible’
An OOPs approach was used for construction in the above animation.
an ant object is a collection of surfaces which have been added together. A unit ant object is shown below.
a colony is a list of ant objects.
a totalAnts object is the list of ants which have been placed at specific locations, then added together to form one object using the get_colony_obj function. The locations are dependent on the frame.
the mobius object is a SurfaceCollection that is geometrically mapped.
The final visualization object is the colony_obj and mobius objects added together to form one object that is then added to the axes.
For the animation, only the colony_obj object needed reconstruction by applying the transform method for each ant copy in the colony list. From a development perspective, the difficultly was determining the Escher_Mobius function which shears the strip around the axial axis of the loop and the corrosponding ant tranformations.
import copy
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.animation as animation
import s3dlib.surface as s3d
import s3dlib.cmap_utilities as cmu
#.. influenced by M.C.Escher - Red Ants, animation
# https://mcescher.com/gallery/mathematical/
# Note: ant objects are rigid and don't bend
# when at the interior of the mobius. MROE
# 0. Define animation control parameters ............................
totalTime, f_domain, numFrames = 4, (0.0,1.0), 60 # time in seconds
frames=np.linspace(*f_domain, numFrames, endpoint=False)
interval = int(1000.0*totalTime/numFrames) # milliseconds
# 1. Define function to examine ....................................
mWidth, mShear = 0.58, 30 # mobius width and shear angle.
def getMesh(Ng,Tu,Tv) :
Su = Ng*Tu # cyclic
Sv = Ng*Tv+1 # bouunded
Uarr = np.linspace(0,1.0,Su+1)
Varr = np.linspace(0,1.0,Sv+1)
y,x = np.meshgrid(Varr,Uarr)
z = np.zeros_like(y)
verts = np.reshape(np.array( [x,y,z] ).T, (-1,3))
fvInd = []
for v in range(0,Sv) :
for u in range(0,Su) :
a = u + v*(Su+1)
d = a + (Su+1)
if not (u%Ng and v%Ng) : fvInd.append( [a,a+1,d+1,d] )
return verts, fvInd
def Escher_Mobius(xyz) :
x,y,z = xyz # xy: [0,1], z=0
U = 2*np.pi*x # domain: [ 0, 2*pi]
V = 2*y - 1 # domain: [-1, 1]
w,phi = mWidth/2, np.pi*mShear/180.0
A = w*np.cos(phi)
B = w*np.sin(phi)
x = w*np.sin(2*U) + A*V*np.sin(U/2)
y = np.sin(U) + B*V*np.sin(U/2)
z = np.cos(U) - w*V*np.cos(U/2)
return x,y,z
def mobius_normals(t) :
w,phi = mWidth/2, np.pi*mShear/180.0
A = w*np.cos(phi)
B = w*np.sin(phi)
U = 2*np.pi*t
V = 0.5 # at mobius midpoint
SS = np.sin(U)*np.sin(U/2) - 0.5*w*V
WC = w*np.cos(U/2)
Nx = -np.cos(U)*WC + B*SS
Ny = 2*w*np.cos(2*U)*WC - A*SS
Nz = 2*w*B*np.cos(2*U)*np.sin(U/2) - A*np.cos(U)*np.sin(U/2)
normal = np.array( [Nx,Ny,Nz])
unitNormal = np.divide(normal, np.linalg.norm(normal,axis=0))
return unitNormal
def loc_to_coor(t) :
def center(t) :
m = 0.5*np.ones(len(t))
z = np.zeros(len(t))
return Escher_Mobius( np.array([t,m,z]) )
h = 0.09
coor = center(2*t) + h*mobius_normals(2*t)
return coor.T
def get_matrix(t) :
t = 2*t
w,phi = mWidth/2, np.pi*mShear/180.0
A = w*np.cos(phi)
B = w*np.sin(phi)
U = 2*np.pi*t
V = 0.5 # at mobius midpoint
Pxu = 2*w*np.cos(2*U) + 0.5*A*V*np.cos(U/2)
Pyu = np.cos(U) + 0.5*B*V*np.cos(U/2)
Pzu = - np.sin(U) + 0.5*w*V*np.sin(U/2)
Pu = np.array( [ Pxu, Pyu, Pzu] )
Pu = np.divide(Pu, np.linalg.norm(Pu,axis=0))
Pxv = A*np.sin(U/2)
Pyv = B*np.sin(U/2)
Pzv = -w*np.cos(U/2)
Pv = np.array( [ Pxv, Pyv, Pzv ] )
Pv = np.divide(Pv, np.linalg.norm(Pv,axis=0))
N = mobius_normals(t)
matrix = np.array( [ Pu, Pv, N ] )
return matrix
def get_ant(length) :
def head_ex(rtp, fX,fY) :
x,y,z = s3d.SphericalSurface.coor_convert(rtp,True)
X = x*np.where(x<0,fX,1)
f,N = 0.5, 3
fact = 1 + f*( ( x/np.amax(x) )**N )
Y = fY*y*np.where(x<0,fact,1)
Z = z*np.where(x<0,fact,1)
return s3d.SphericalSurface.coor_convert([X,Y,Z],False)
def stretch(rtp ,f) :
x,y,z = s3d.SphericalSurface.coor_convert(rtp,True)
return s3d.SphericalSurface.coor_convert([f*x,y,z],False)
def abdomen_ex(rtp, fZ) :
x,y,z = s3d.SphericalSurface.coor_convert(rtp,True)
f,N = 0.5, 3
fact = 1 - f*( ( z/np.amax(z) )**N )
Z = z*np.where(z>0,fZ,1)
Y = y*np.where(z>0,fact,1)
X = x*np.where(z>0,fact,1)
return s3d.SphericalSurface.coor_convert([X,Y,Z],False)
# ..............................................................
rez,antcolor,antcolor2 = 2, [0.556, 0.266, 0.206], [0.898, 0.686, 0.561]
bw = cmu.stitch_color('bisque',antcolor, 'k',antcolor,'bisque',bndry=[.07,0.15,0.85,.93])
ag = cmu.rgb_cmap_gradient(antcolor,antcolor2)
abdomen = s3d.SphericalSurface.grid(rez*5,rez*10).domain(22)
abdomen.map_geom_from_op(lambda c: abdomen_ex(c,1.8) )
abdomen.transform( s3d.eulerRot(0,90+20,useXconv=False), translate=[40,0,11] )
bridge = s3d.CylindricalSurface(rez-1).domain(4,20)
bridge.transform( s3d.eulerRot(0,90,useXconv=False), translate=[25,0,0] )
thorax_r = s3d.SphericalSurface(rez).domain(8.0)
thorax_r.map_geom_from_op( lambda c: stretch(c,1.900) )
thorax_f = s3d.SphericalSurface(rez).domain(10.5)
thorax_f.map_geom_from_op( lambda c: stretch(c,1.600))
thorax_f.transform( s3d.eulerRot(0,40,useXconv=False), translate=[-22,0,3] )
body = thorax_f + thorax_r + bridge + abdomen
body.map_cmap_from_normals(ag,direction=[0,0,1])
head = s3d.SphericalSurface(rez,color=antcolor).domain(13.5)
head.map_geom_from_op( lambda c: head_ex(c, 1.99, 1.4) )
head.map_cmap_from_normals(ag,direction=[0,0,1])
eye_r = s3d.SphericalSurface(rez+1,color='k').domain(7).transform(translate=[0,-16,0])
eye_r.map_cmap_from_op(lambda c: s3d.SphericalSurface.coor_convert(c,True)[1]**4, bw)
eye_l = s3d.SphericalSurface(rez+1,color='k').domain(7).transform(translate=[0, 16,0])
eye_l.map_cmap_from_op(lambda c: s3d.SphericalSurface.coor_convert(c,True)[1]**4, bw)
head += eye_r + eye_l
head.transform( s3d.eulerRot(0,-51,useXconv=False), translate=[-48,0,8] )
ant = head + body
scale = 1/(ant.bounds['xlim'][1] - ant.bounds['xlim'][0] )
ant.transform(scale=length*scale)
return ant
def get_colony_object(frm,arrayOfAnts) :
numbAnts = len(arrayOfAnts)
colonyCopy = copy.deepcopy(arrayOfAnts)
initial_loc = np.linspace(0,1,numbAnts,endpoint=False)
current_loc = initial_loc - frm/numbAnts
current_coor = loc_to_coor(current_loc)
totalAnts = None
for i, ant in enumerate(colonyCopy) :
matrix = get_matrix(current_loc[i])
ant.transform(matrix,translate=current_coor[i])
if totalAnts is None : totalAnts = ant
else: totalAnts += ant
totalAnts.shade(direction=[-1,1,1]).hilite(0.8,focus=2,direction=[-1,1,1])
return totalAnts
# 2. Setup and map surface .........................................
gridcolor = [0.788,0.729,0.647]
Ng,Tu,Tv = 5, 50, 4 # grid layout parameters
xyz,fvInd = getMesh(Ng,Tu,Tv)
mobius = s3d.Surface3DCollection(xyz, fvInd, color=gridcolor).triangulate()
mobius.map_geom_from_op(Escher_Mobius)
mobius.shade(0.3,direction=[1,-1,1], isAbs=True)
numbAnts, antlength = 9, 0.58
colony = [ get_ant(antlength) for i in range(numbAnts) ]
frame = 0.6
# move colony to the frame location.
totalAnts = get_colony_object(frame,colony)
redAnts = mobius + totalAnts
# 3. Construct figure, add surface, and plot .......................
info,textColor = 'S3Dlib.org', [0.502, 0.200, 0.278, 0.3]
minmax = (-0.8,0.8) # axis limits
fig = plt.figure(figsize=(6,7), constrained_layout=True, facecolor='k')
text = fig.text(0.97, 0.14, info, color=textColor, ha='right',
rotation=90, fontsize=70, fontweight='bold' )
ax = plt.axes(projection='3d', facecolor='k', proj_type='ortho')
ax.set(xlim=minmax, ylim=minmax, zlim=minmax)
ax.set_axis_off()
ax.view_init(89.9,-89.9)
ax.add_collection3d(redAnts)
#plt.show()
# 4. Animation ......................................................
def update_fig(frame):
global redAnts
global colony
redAnts.remove()
totalAnts = get_colony_object(frame,colony)
redAnts = mobius + totalAnts
ax.add_collection3d(redAnts)
return redAnts,
anim = FuncAnimation(fig, update_fig, frames, interval=interval, repeat=True)
anim.save('redants.html',writer='html')
msg = "saved {} frames, values: [{:.3f} to {:.3f}] @ {} milliseconds/frane"
print(msg.format(numFrames,np.min(frames),np.max(frames),interval))