"""Refine the existing Blender slab, preserving a shared regular height lattice."""
import bpy, math, json
import numpy as np
from pathlib import Path
from mathutils import Vector
from mathutils.bvhtree import BVHTree
ROOT=Path('D:/project/sunrise')
s=bpy.data.scenes['Summit slab and marker'];bpy.context.window.scene=s
slab=bpy.data.objects['SLAB_SURFACE']
if slab.get('jointStudyVersion'):raise RuntimeError('Study already refined; edit it in place.')
backup=slab.copy();backup.data=slab.data.copy();backup.name='Slab before joint refinement';s.collection.objects.link(backup);backup.hide_render=True;backup.hide_set(True)
bvh=BVHTree.FromObject(slab,bpy.context.evaluated_depsgraph_get())
outline=bpy.data.objects['WALKABLE_OUTLINE'];poly=[(v.co.x,v.co.y) for v in outline.data.vertices]
joints=[ [(-26,7),(-18,10),(-9,12),(0,11),(12,15)], [(-26,22),(-17,25),(-8,26),(1,29),(14,28)], [(-24,39),(-14,40),(-4,43),(6,42),(14,46)], [(-20,54),(-10,51),(0,54),(11,53)], [(-18,-7),(-15,4),(-11,12),(-14,23),(-9,33),(-11,43),(-7,57)], [(3,-9),(1,3),(3,11),(-1,21),(2,29),(0,40),(5,53)] ]
def distances(x,y,lines):
    out=np.full(np.broadcast(x,y).shape,1e5)
    for line in lines:
        for (ax,ay),(bx,by) in zip(line,line[1:]):
            dx=bx-ax;dy=by-ay;t=np.clip(((x-ax)*dx+(y-ay)*dy)/(dx*dx+dy*dy),0,1)
            out=np.minimum(out,np.hypot(x-ax-t*dx,y-ay-t*dy))
    return out

def edge_distance(x,y):
    d=distances(x,y,[poly+[poly[0]]]);hit=np.zeros(np.broadcast(x,y).shape,dtype=bool)
    for (ax,ay),(bx,by) in zip(poly,poly[1:]+poly[:1]):
        if ay!=by:hit^=((ay>y)!=(by>y))&(x<(bx-ax)*(y-ay)/(by-ay)+ax)
    return np.where(hit,-d,d)

def sculpt(x,y):
    d=edge_distance(x,y);fade=1-np.clip((d+1)/5,0,1);fade=fade*fade*(3-2*fade)
    cracks=distances(x,y,joints)
    roll=.72*np.sin(x*.19+y*.055)*np.cos(y*.13)+.36*np.sin(y*.33+x*.09)
    shoulders=.42*np.tanh((y-18-.18*x)/.85)+.32*np.tanh((y-38+.13*x)/.7)
    cut=.50*np.exp(-(cracks/.38)**2)+.15*np.exp(-(cracks/1.2)**2)
    lip=-.7*np.exp(-((d+.25)/1.1)**2)
    return fade*(roll+shoulders-cut+lip)
# Ten subdivisions per original 5m interval. Retain every original XY coordinate.
xs=sorted(set(round(v.co.x,6) for v in slab.data.vertices));ys=sorted(set(round(v.co.y,6) for v in slab.data.vertices))
xs=[a+(b-a)*j/10 for a,b in zip(xs,xs[1:]) for j in range(10)]+[xs[-1]]
ys=[a+(b-a)*j/10 for a,b in zip(ys,ys[1:]) for j in range(10)]+[ys[-1]]
verts=[]
for y in ys:
    for x in xs:
        hit=bvh.ray_cast(Vector((x,y,100)),Vector((0,0,-1)))[0]
        if hit is None:hit=bvh.find_nearest(Vector((x,y,0)))[0]
        verts.append((x,y,hit.z+float(sculpt(np.asarray(x),np.asarray(y)))))
nx=len(xs);faces=[]
for j in range(len(ys)-1):
    for i in range(nx-1):
        a=j*nx+i;faces.extend([(a,a+1,a+nx),(a+1,a+nx+1,a+nx)])
mesh=bpy.data.meshes.new('Sculpted jointed granite lattice');mesh.from_pydata(verts,[],faces);mesh.update();slab.data=mesh
# Unique atlas: warm granite grains, large weathering patches and matching fissures.
N=1536
X,Y=np.meshgrid(np.linspace(xs[0],xs[-1],N),np.linspace(ys[0],ys[-1],N))
D=distances(X,Y,joints)
rng=np.random.default_rng(720)
weather=.035*np.sin(X*.41+np.sin(Y*.17))+.025*np.cos(Y*.57-X*.13)
grain=rng.normal(0,.018,(N,N));fleck=(rng.random((N,N))>.96)*-.055
crack=.22*np.exp(-(D/.16)**2)+.065*np.exp(-(D/.7)**2)
rgb=np.stack([.64+weather+grain+fleck-crack,.59+weather+grain+fleck-crack,.49+weather+grain+fleck-crack],axis=-1)
rgba=np.concatenate([np.clip(rgb,0,1),np.ones((N,N,1))],axis=2).astype(np.float32)
im=bpy.data.images.new('Summit unique weathered granite',width=N,height=N,alpha=True);im.pixels.foreach_set(rgba.ravel());im.update()
im.filepath_raw=str(ROOT/'public/terrain/materials/plateau-color.png')
im.file_format=next(i.identifier for i in im.bl_rna.properties['file_format'].enum_items if i.identifier=='PNG');im.save()
mat=bpy.data.materials.new('Summit unique granite atlas');mat.use_nodes=True
bsdf=next(n for n in mat.node_tree.nodes if n.type=='BSDF_PRINCIPLED');bsdf.inputs['Roughness'].default_value=.92
tex=mat.node_tree.nodes.new('ShaderNodeTexImage');tex.image=im;mat.node_tree.links.new(tex.outputs['Color'],bsdf.inputs['Base Color']);mesh.materials.append(mat)
uv=mesh.uv_layers.new(name='Unique plateau atlas')
for p in mesh.polygons:
    p.use_smooth=True
    for li in p.loop_indices:
        v=mesh.vertices[mesh.loops[li].vertex_index].co;uv.data[li].uv=((v.x-xs[0])/(xs[-1]-xs[0]),(v.y-ys[0])/(ys[-1]-ys[0]))
slab['jointStudyVersion']=1
slab['atlasBounds']=[xs[0],-ys[-1],xs[-1],-ys[0]]
marker=bpy.data.objects['summit-marker'];marker.location.z+=float(sculpt(np.asarray(marker.location.x),np.asarray(marker.location.y)))
print('Sculpted',len(verts),'vertices;',len(faces),'triangles; six intersecting fracture traces and weathered edge lip')