# Creator 1 mesh iteration: "Threads Give Way" - a plain-weave cloth (warp+weft threads) pinned at the frame,
# a sphere is pushed through the middle. Position-based dynamics: segment-length constraints along threads,
# crossing (shear) links between warp and weft, bending via skip-links, approximate self-contact via KD-tree,
# sphere collision opens a hole whose rim is projected. Colors: Flexoki.
import numpy as np
from scipy.spatial import cKDTree
from PIL import Image, ImageDraw
rng=np.random.default_rng(11)
N=44                     # threads per direction
S=1.0/(N-1)
u,v=np.meshgrid(np.linspace(-.5,.5,N),np.linspace(-.5,.5,N),indexing='ij')
P=np.stack([u,v,np.zeros_like(u)],-1).reshape(-1,3); P[:,2]+=rng.normal(0,1e-4,len(P))
ID=lambda i,j:i*N+j
rest=P.copy(); Pr=P.copy()
cons=[]
for i in range(N):
    for j in range(N):
        if i+1<N: cons.append((ID(i,j),ID(i+1,j),S))
        if j+1<N: cons.append((ID(i,j),ID(i,j+1),S))
        if i+1<N and j+1<N:
            cons.append((ID(i,j),ID(i+1,j+1),S*2**.5)); cons.append((ID(i+1,j),ID(i,j+1),S*2**.5))
        if i+2<N: cons.append((ID(i,j),ID(i+2,j),2*S*1.0))
        if j+2<N: cons.append((ID(i,j),ID(i,j+2),2*S*1.0))
C=np.array([(a,b) for a,b,_ in cons]); L=np.array([l for *_,l in cons])
edge=(np.abs(rest[:,0])>.499)|(np.abs(rest[:,1])>.499)
sph_c=np.array([0.04,-0.03,0.0]); R=0.17
V=np.zeros_like(P)
cnt=np.bincount(C.ravel(),minlength=len(P))[:,None].astype(float)
for step in range(260):
    t=min(1,step/180)
    sph_c[2]=-0.30+0.62*t            # sphere travels through the plane (towards viewer)
    Pn=P+V*0.5+np.array([0,0,0])
    Pn[:,2]+=0  # no gravity
    # sphere collision (position projection)
    for it in range(14):
        d=Pn[C[:,1]]-Pn[C[:,0]]; l=np.linalg.norm(d,axis=1)+1e-9
        k=((l-L)/l)[:,None]*d*0.5
        # compression/shear stiff, stretch softer
        D=np.zeros_like(Pn); np.add.at(D,C[:,0],k); np.add.at(D,C[:,1],-k); Pn+=D*1.6/cnt
        r=Pn-sph_c; rd=np.linalg.norm(r,axis=1)
        m=rd<R*1.03
        Pn[m]=sph_c+r[m]/rd[m,None]*R*1.03
        tree=cKDTree(Pn); pr=tree.query_pairs(S*0.55,output_type='ndarray')
        if len(pr):
            d=Pn[pr[:,1]]-Pn[pr[:,0]]; l=np.linalg.norm(d,axis=1)+1e-9
            ok=(np.abs(pr[:,0]-pr[:,1])>2)&(np.abs(pr[:,0]//N-pr[:,1]//N)+np.abs(pr[:,0]%N-pr[:,1]%N)>2)
            pr=pr[ok];d=d[ok];l=l[ok]
            k=((S*.55-l)/l)[:,None]*d*0.5
            np.add.at(Pn,pr[:,0],-k); np.add.at(Pn,pr[:,1],k)
        Pn[edge]=rest[edge]
    V=np.clip(Pn-P,-.02,.02); P=Pn
# projected opening: the sphere's silhouette rim in view space
ang=np.radians(24); ca,sa=np.cos(ang),np.sin(ang)
def proj(p):
    x=p[...,0]; y=p[...,1]*ca-p[...,2]*sa; z=p[...,1]*sa+p[...,2]*ca
    return x,y,z
x,y,z=proj(P)
W=1200; M=90; sc=W-2*M
def T(x,y): return (M+(x+.5)*sc, M+(y+.5)*sc)
BG=(16,15,15); cols=[(209,77,65),(218,112,44),(208,162,21),(135,154,57),(58,169,159),(67,133,190),(139,126,200),(206,93,151)]
img=Image.new("RGB",(W,W),BG); d=ImageDraw.Draw(img)
# opening: dark ellipse of sphere projection + bright rim
cx,cy,_=proj(sph_c); rr=R*sc
c0=T(cx,cy)
d.ellipse([c0[0]-rr*1.0,c0[1]-rr*1.0,c0[0]+rr*1.0,c0[1]+rr*1.0],fill=(8,8,8))
segs=[]
for i in range(N):
    for j in range(N):
        a=ID(i,j)
        if i+1<N: segs.append((a,ID(i+1,j),'w',i))
        if j+1<N: segs.append((a,ID(i,j+1),'f',j))
def draw_seg(a,b,kind,k,thick):
    zz=(z[a]+z[b])/2
    # plain weave: threads alternate over/under by parity, modulated by depth
    par=((a//N+a%N)%2==0)
    zz+= (0.004 if (par==(kind=='w')) else -0.004)
    return zz
items=[]
for a,b,kind,k in segs:
    items.append((draw_seg(a,b,kind,k,0),a,b,kind,k))
items.sort(key=lambda q:q[0])
for zz,a,b,kind,k in items:
    stretch=np.linalg.norm(P[a]-P[b])/S
    # color by distance from the sphere axis (how far the thread was pushed), bands
    dist=np.linalg.norm(((P[a]+P[b])/2)[:2]-sph_c[:2])
    ci=int(np.clip(dist/0.62*len(cols),0,len(cols)-1))
    c=cols[ci] if kind=='w' else cols[(ci+4)%len(cols)]
    f=np.clip(0.55+ (zz+0.15)*1.6,0.35,1.0)
    c=tuple(int(BG[q]+(c[q]-BG[q])*f) for q in range(3))
    d.line([T(x[a],y[a]),T(x[b],y[b])],fill=c,width=3)
# opening rim
for w_,col in [(5,(206,205,195))]:
    d.ellipse([c0[0]-rr,c0[1]-rr,c0[0]+rr,c0[1]+rr],outline=col,width=2)

FLEXOKI="fffcf0 f2f0e5 e6e4d9 cecdc3 b7b5ac 9f9d96 878580 6f6e69 575653 403e3c 343331 282726 1c1b1a 100f0f d14d41 da702c d0a215 879a39 3aa99f 4385be 8b7ec8 ce5d97 af3029 bc5215 ad8301 66800b 24837b 205ea6 5e409d a02f6f".split()
_pal=np.array([[int(h[i:i+2],16) for i in (0,2,4)] for h in FLEXOKI],dtype=np.int32)
_a=np.asarray(img.convert('RGB')).astype(np.int32)
_flat=_a.reshape(-1,3)
_out=np.empty(len(_flat),dtype=np.int32)
for _s in range(0,len(_flat),200000):
    _d=((_flat[_s:_s+200000,None,:]-_pal[None])**2).sum(-1)
    _out[_s:_s+200000]=_d.argmin(1)
img=Image.fromarray(_pal[_out].reshape(_a.shape).astype(np.uint8))
img.save("render.png")