```python
import bpy, bmesh, math, random
from mathutils import Vector, Matrix, Euler
from mathutils import noise as mnoise

random.seed(11)
TAU = math.pi * 2.0

# ----------------------------------------------------------------------------- clean
for ob in list(bpy.data.objects):
    bpy.data.objects.remove(ob, do_unlink=True)
for me in list(bpy.data.meshes):
    bpy.data.meshes.remove(me)
for mt in list(bpy.data.materials):
    bpy.data.materials.remove(mt)

scene = bpy.context.scene

# ----------------------------------------------------------------------------- helpers
def link(ob):
    scene.collection.objects.link(ob)
    return ob

def obj_from_bm(name, bm, mat, smooth=False):
    me = bpy.data.meshes.new(name)
    bm.to_mesh(me)
    bm.free()
    if mat:
        me.materials.append(mat)
    if smooth:
        for p in me.polygons:
            p.use_smooth = True
    ob = bpy.data.objects.new(name, me)
    return link(ob)

def add_box(bm, M, size):
    bmesh.ops.create_cube(bm, size=1.0,
        matrix=M @ Matrix.Diagonal(Vector((size[0], size[1], size[2], 1.0))))

def add_tube(bm, pts, radius, sides=6, radii=None):
    n = len(pts)
    if n < 2:
        return
    rings = []
    for i, p in enumerate(pts):
        if i < n - 1:
            d = pts[i + 1] - pts[i]
        else:
            d = pts[i] - pts[i - 1]
        if d.length < 1e-6:
            d = Vector((0, 0, 1))
        d = d.normalized()
        up = Vector((0, 0, 1))
        if abs(d.dot(up)) > 0.95:
            up = Vector((1, 0, 0))
        ax = d.cross(up).normalized()
        ay = d.cross(ax).normalized()
        r = radius if radii is None else max(0.008, radii[i])
        ring = []
        for k in range(sides):
            a = TAU * k / sides
            ring.append(bm.verts.new(p + ax * (math.cos(a) * r) + ay * (math.sin(a) * r)))
        rings.append(ring)
    for i in range(n - 1):
        for k in range(sides):
            k2 = (k + 1) % sides
            try:
                bm.faces.new([rings[i][k], rings[i][k2], rings[i + 1][k2], rings[i + 1][k]])
            except Exception:
                pass
    try:
        bm.faces.new(rings[0][::-1])
    except Exception:
        pass
    try:
        bm.faces.new(rings[-1])
    except Exception:
        pass

LEAF_PTS = [Vector((0, 0, 0)), Vector((0.34, 0.42, 0.06)),
            Vector((0, 1.0, 0.0)), Vector((-0.34, 0.42, 0.06))]

def rand_rot(rnd):
    return Euler((rnd.uniform(0, TAU), rnd.uniform(0, TAU), rnd.uniform(0, TAU)), 'XYZ').to_matrix().to_4x4()

def add_leaf(bm, center, size, rot):
    M = Matrix.Translation(center) @ rot @ Matrix.Diagonal(Vector((size, size, size, 1.0)))
    try:
        bm.faces.new([bm.verts.new(M @ p) for p in LEAF_PTS])
    except Exception:
        pass

def leaf_cluster(bm, center, radius, count, size, rnd, flatten=1.0):
    for _ in range(count):
        d = Vector((rnd.gauss(0, 1), rnd.gauss(0, 1), rnd.gauss(0, 1) * flatten))
        if d.length < 1e-5:
            d = Vector((0, 0, 1))
        p = Vector(center) + d.normalized() * (radius * rnd.uniform(0.35, 1.0))
        add_leaf(bm, p, size * rnd.uniform(0.65, 1.35), rand_rot(rnd))

def add_blade(bm, base, height, ang, bend, width):
    dirv = Vector((math.cos(ang), math.sin(ang), 0))
    perp = Vector((-dirv.y, dirv.x, 0))
    rows = []
    N = 3
    for i in range(N + 1):
        t = i / N
        z = height * (t - 0.12 * t * t)
        p = Vector(base) + dirv * (bend * height * t * t) + Vector((0, 0, z))
        hw = width * (1 - 0.9 * t)
        rows.append((p - perp * hw, p + perp * hw))
    for i in range(N):
        a, b = rows[i]
        c, dd = rows[i + 1]
        try:
            bm.faces.new([bm.verts.new(a), bm.verts.new(b), bm.verts.new(dd), bm.verts.new(c)])
        except Exception:
            pass

def add_corrugated(bm, M, lx, ly, amp=0.035, period=0.30, nv=4, sag=0.0, y0=0.0):
    nu = max(6, int(lx / period * 3))
    grid = []
    for i in range(nu + 1):
        col = []
        x = -lx * 0.5 + lx * i / nu
        cz = amp * math.sin(TAU * x / period)
        for j in range(nv + 1):
            y = y0 + (ly - y0) * j / nv
            tv = (y - y0) / max(1e-6, (ly - y0))
            z = cz - sag * math.sin(math.pi * tv)
            col.append(bm.verts.new(M @ Vector((x, y, z))))
        grid.append(col)
    for i in range(nu):
        for j in range(nv):
            try:
                bm.faces.new([grid[i][j], grid[i + 1][j], grid[i + 1][j + 1], grid[i][j + 1]])
            except Exception:
                pass

# ----------------------------------------------------------------------------- materials
def base_mat(name):
    m = bpy.data.materials.new(name)
    m.use_nodes = True
    nt = m.node_tree
    for n in list(nt.nodes):
        nt.nodes.remove(n)
    out = nt.nodes.new('ShaderNodeOutputMaterial'); out.location = (700, 0)
    b = nt.nodes.new('ShaderNodeBsdfPrincipled'); b.location = (380, 0)
    nt.links.new(b.outputs[0], out.inputs[0])
    return m, nt, b

def sset(node, key, val):
    s = node.inputs.get(key)
    if s is not None:
        try:
            s.default_value = val
        except Exception:
            pass

def mix_node(nt, loc=(0, 0)):
    n = nt.nodes.new('ShaderNodeMix')
    n.data_type = 'RGBA'
    n.blend_type = 'MIX'
    n.location = loc
    return n

def mix_set(n, a=None, b=None, fac=None):
    if a is not None:
        n.inputs[6].default_value = a
    if b is not None:
        n.inputs[7].default_value = b
    if fac is not None:
        n.inputs[0].default_value = fac

def mix_out(n):
    return n.outputs[2]

def noise_node(nt, loc, scale, detail=6.0, rough=0.55, vector=None, links_from=None):
    n = nt.nodes.new('ShaderNodeTexNoise')
    n.location = loc
    sset(n, 'Scale', scale); sset(n, 'Detail', detail); sset(n, 'Roughness', rough)
    if links_from is not None:
        nt.links.new(links_from, n.inputs['Vector'])
    return n

def ramp(nt, loc, stops):
    n = nt.nodes.new('ShaderNodeValToRGB')
    n.location = loc
    cr = n.color_ramp
    cr.elements[0].position = stops[0][0]
    cr.elements[0].color = stops[0][1]
    cr.elements[1].position = stops[-1][0]
    cr.elements[1].color = stops[-1][1]
    for pos, col in stops[1:-1]:
        e = cr.elements.new(pos)
        e.color = col
    return n

def mat_rust():
    m, nt, b = base_mat("RustedMetalRoof")
    tc = nt.nodes.new('ShaderNodeTexCoord'); tc.location = (-900, 0)
    n1 = noise_node(nt, (-700, 100), 2.6, 9.0, 0.62, links_from=tc.outputs['Object'])
    n2 = noise_node(nt, (-700, -250), 14.0, 7.0, 0.55, links_from=tc.outputs['Object'])
    col = ramp(nt, (-450, 150), [
        (0.25, (0.040, 0.024, 0.017, 1)),
        (0.40, (0.105, 0.048, 0.024, 1)),
        (0.52, (0.230, 0.090, 0.032, 1)),
        (0.63, (0.330, 0.140, 0.050, 1)),
        (0.78, (0.150, 0.135, 0.120, 1)),
        (0.92, (0.230, 0.225, 0.215, 1)),
    ])
    nt.links.new(n1.outputs['Fac'], col.inputs['Fac'])
    # streak modulation
    mx = mix_node(nt, (-180, 120))
    mx.blend_type = 'MULTIPLY'
    mix_set(mx, fac=0.35)
    nt.links.new(col.outputs['Color'], mx.inputs[6])
    nt.links.new(n2.outputs['Color'], mx.inputs[7])
    nt.links.new(mix_out(mx), b.inputs['Base Color'])
    rg = ramp(nt, (-450, -80), [(0.35, (0.98, 0.98, 0.98, 1)), (0.95, (0.42, 0.42, 0.42, 1))])
    nt.links.new(n1.outputs['Fac'], rg.inputs['Fac'])
    nt.links.new(rg.outputs['Color'], b.inputs['Roughness'])
    mg = ramp(nt, (-450, -320), [(0.72, (0, 0, 0, 1)), (0.95, (0.7, 0.7, 0.7, 1))])
    nt.links.new(n1.outputs['Fac'], mg.inputs['Fac'])
    nt.links.new(mg.outputs['Color'], b.inputs['Metallic'])
    nb = noise_node(nt, (-700, -520), 60.0, 8.0, 0.6, links_from=tc.outputs['Object'])
    bump = nt.nodes.new('ShaderNodeBump'); bump.location = (100, -350)
    sset(bump, 'Strength', 0.45); sset(bump, 'Distance', 0.03)
    nt.links.new(nb.outputs['Fac'], bump.inputs['Height'])
    nt.links.new(bump.outputs['Normal'], b.inputs['Normal'])
    return m

def mat_wall():
    m, nt, b = base_mat("WeatheredWall")
    tc = nt.nodes.new('ShaderNodeTexCoord'); tc.location = (-1100, 0)
    mp = nt.nodes.new('ShaderNodeMapping'); mp.location = (-950, -200)
    mp.inputs['Scale'].default_value = (3.0, 3.0, 0.5)
    nt.links.new(tc.outputs['Object'], mp.inputs['Vector'])
    stain = noise_node(nt, (-760, -200), 5.0, 9.0, 0.65, links_from=mp.outputs['Vector'])
    plaster = mix_node(nt, (-520, -120))
    mix_set(plaster, a=(0.330, 0.300, 0.262, 1), b=(0.130, 0.115, 0.100, 1))
    nt.links.new(stain.outputs['Fac'], plaster.inputs[0])
    # moss rising from the ground
    sep = nt.nodes.new('ShaderNodeSeparateXYZ'); sep.location = (-950, 250)
    nt.links.new(tc.outputs['Generated'], sep.inputs[0])
    mn = noise_node(nt, (-950, 480), 6.0, 6.0, 0.6, links_from=tc.outputs['Object'])
    mul = nt.nodes.new('ShaderNodeMath'); mul.location = (-760, 480)
    mul.operation = 'MULTIPLY'; mul.inputs[1].default_value = 0.45
    nt.links.new(mn.outputs['Fac'], mul.inputs[0])
    add = nt.nodes.new('ShaderNodeMath'); add.location = (-600, 350)
    add.operation = 'ADD'
    nt.links.new(sep.outputs['Z'], add.inputs[0])
    nt.links.new(mul.outputs[0], add.inputs[1])
    mramp = ramp(nt, (-420, 350), [(0.10, (1, 1, 1, 1)), (0.62, (0, 0, 0, 1))])
    nt.links.new(add.outputs[0], mramp.inputs['Fac'])
    mossmix = mix_node(nt, (-160, 100))
    mix_set(mossmix, b=(0.055, 0.095, 0.032, 1))
    nt.links.new(mramp.outputs['Color'], mossmix.inputs[0])
    nt.links.new(mix_out(plaster), mossmix.inputs[6])
    nt.links.new(mix_out(mossmix), b.inputs['Base Color'])
    sset(b, 'Roughness', 0.92)
    nb = noise_node(nt, (-760, -520), 22.0, 8.0, 0.6, links_from=tc.outputs['Object'])
    bump = nt.nodes.new('ShaderNodeBump'); bump.location = (100, -350)
    sset(bump, 'Strength', 0.45); sset(bump, 'Distance', 0.05)
    nt.links.new(nb.outputs['Fac'], bump.inputs['Height'])
    nt.links.new(bump.outputs['Normal'], b.inputs['Normal'])
    return m

def mat_wood():
    m, nt, b = base_mat("RottenWood")
    tc = nt.nodes.new('ShaderNodeTexCoord'); tc.location = (-800, 0)
    mp = nt.nodes.new('ShaderNodeMapping'); mp.location = (-650, 0)
    mp.inputs['Scale'].default_value = (14.0, 1.2, 1.2)
    nt.links.new(tc.outputs['Object'], mp.inputs['Vector'])
    n = noise_node(nt, (-470, 0), 6.0, 9.0, 0.7, links_from=mp.outputs['Vector'])
    c = ramp(nt, (-260, 0), [(0.30, (0.038, 0.026, 0.018, 1)),
                             (0.55, (0.085, 0.058, 0.036, 1)),
                             (0.85, (0.140, 0.105, 0.070, 1))])
    nt.links.new(n.outputs['Fac'], c.inputs['Fac'])
    nt.links.new(c.outputs['Color'], b.inputs['Base Color'])
    sset(b, 'Roughness', 0.9)
    bump = nt.nodes.new('ShaderNodeBump'); bump.location = (120, -300)
    sset(bump, 'Strength', 0.4)
    nt.links.new(n.outputs['Fac'], bump.inputs['Height'])
    nt.links.new(bump.outputs['Normal'], b.inputs['Normal'])
    return m

def mat_leaf():
    m, nt, b = base_mat("Foliage")
    tc = nt.nodes.new('ShaderNodeTexCoord'); tc.location = (-700, 0)
    n = noise_node(nt, (-520, 0), 0.7, 5.0, 0.55, links_from=tc.outputs['Object'])
    c = ramp(nt, (-300, 0), [(0.30, (0.020, 0.055, 0.012, 1)),
                             (0.48, (0.055, 0.120, 0.026, 1)),
                             (0.66, (0.100, 0.180, 0.040, 1)),
                             (0.86, (0.170, 0.170, 0.045, 1))])
    nt.links.new(n.outputs['Fac'], c.inputs['Fac'])
    nt.links.new(c.outputs['Color'], b.inputs['Base Color'])
    sset(b, 'Roughness', 0.55)
    sset(b, 'Subsurface Weight', 0.12)
    sset(b, 'Subsurface Radius', (0.03, 0.06, 0.02))
    return m

def mat_grass():
    m, nt, b = base_mat("Grass")
    tc = nt.nodes.new('ShaderNodeTexCoord'); tc.location = (-700, 0)
    sep = nt.nodes.new('ShaderNodeSeparateXYZ'); sep.location = (-520, 120)
    nt.links.new(tc.outputs['Generated'], sep.inputs[0])
    c = ramp(nt, (-320, 120), [(0.05, (0.030, 0.070, 0.014, 1)),
                               (0.55, (0.090, 0.150, 0.032, 1)),
                               (1.00, (0.175, 0.170, 0.055, 1))])
    nt.links.new(sep.outputs['Z'], c.inputs['Fac'])
    nt.links.new(c.outputs['Color'], b.inputs['Base Color'])
    sset(b, 'Roughness', 0.6)
    sset(b, 'Subsurface Weight', 0.1)
    return m

def mat_ground():
    m, nt, b = base_mat("MossyGround")
    tc = nt.nodes.new('ShaderNodeTexCoord'); tc.location = (-1200, 0)
    sep = nt.nodes.new('ShaderNodeSeparateXYZ'); sep.location = (-1030, 200)
    nt.links.new(tc.outputs['Object'], sep.inputs[0])
    ab = nt.nodes.new('ShaderNodeMath'); ab.location = (-870, 200)
    ab.operation = 'ABSOLUTE'
    nt.links.new(sep.outputs['X'], ab.inputs[0])
    sc = nt.nodes.new('ShaderNodeMath'); sc.location = (-720, 200)
    sc.operation = 'DIVIDE'; sc.inputs[1].default_value = 7.0
    nt.links.new(ab.outputs[0], sc.inputs[0])
    pn = noise_node(nt, (-1030, -60), 1.4, 8.0, 0.6, links_from=tc.outputs['Object'])
    sub = nt.nodes.new('ShaderNodeMath'); sub.location = (-870, -60)
    sub.operation = 'SUBTRACT'; sub.inputs[1].default_value = 0.5
    nt.links.new(pn.outputs['Fac'], sub.inputs[0])
    mu = nt.nodes.new('ShaderNodeMath'); mu.location = (-720, -60)
    mu.operation = 'MULTIPLY'; mu.inputs[1].default_value = 0.30
    nt.links.new(sub.outputs[0], mu.inputs[0])
    ad = nt.nodes.new('ShaderNodeMath'); ad.location = (-560, 100)
    ad.operation = 'ADD'
    nt.links.new(sc.outputs[0], ad.inputs[0])
    nt.links.new(mu.outputs[0], ad.inputs[1])
    r = ramp(nt, (-390, 100), [(0.20, (0, 0, 0, 1)), (0.50, (1, 1, 1, 1))])
    nt.links.new(ad.outputs[0], r.inputs['Fac'])
    mx = mix_node(nt, (-160, 60))
    mix_set(mx, a=(0.115, 0.088, 0.062, 1), b=(0.055, 0.090, 0.028, 1))
    nt.links.new(r.outputs['Color'], mx.inputs[0])
    n2 = noise_node(nt, (-390, -320), 3.5, 8.0, 0.65, links_from=tc.outputs['Object'])
    mx2 = mix_node(nt, (60, 0))
    mix_set(mx2, b=(0.095, 0.125, 0.045, 1), fac=0.5)
    nt.links.new(mix_out(mx), mx2.inputs[6])
    nt.links.new(n2.outputs['Fac'], mx2.inputs[0])
    nt.links.new(mix_out(mx2), b.inputs['Base Color'])
    sset(b, 'Roughness', 0.93)
    nb = noise_node(nt, (-390, -560), 26.0, 8.0, 0.6, links_from=tc.outputs['Object'])
    bump = nt.nodes.new('ShaderNodeBump'); bump.location = (200, -350)
    sset(bump, 'Strength', 0.55); sset(bump, 'Distance', 0.08)
    nt.links.new(nb.outputs['Fac'], bump.inputs['Height'])
    nt.links.new(bump.outputs['Normal'], b.inputs['Normal'])
    return m

def mat_stone():
    m, nt, b = base_mat("Stone")
    tc = nt.nodes.new('ShaderNodeTexCoord'); tc.location = (-700, 0)
    n = noise_node(nt, (-520, 0), 8.0, 8.0, 0.65, links_from=tc.outputs['Object'])
    c = ramp(nt, (-300, 0), [(0.30, (0.075, 0.072, 0.065, 1)),
                             (0.60, (0.140, 0.135, 0.122, 1)),
                             (0.85, (0.070, 0.095, 0.045, 1))])
    nt.links.new(n.outputs['Fac'], c.inputs['Fac'])
    nt.links.new(c.outputs['Color'], b.inputs['Base Color'])
    sset(b, 'Roughness', 0.95)
    bump = nt.nodes.new('ShaderNodeBump'); bump.location = (120, -300)
    sset(bump, 'Strength', 0.6); sset(bump, 'Distance', 0.06)
    nt.links.new(n.outputs['Fac'], bump.inputs['Height'])
    nt.links.new(bump.outputs['Normal'], b.inputs['Normal'])
    return m

def mat_brick():
    m, nt, b = base_mat("OldBrick")
    tc = nt.nodes.new('ShaderNodeTexCoord'); tc.location = (-800, 0)
    br = nt.nodes.new('ShaderNodeTexBrick'); br.location = (-600, 0)
    br.inputs['Scale'].default_value = 9.0
    br.inputs['Color1'].default_value = (0.140, 0.062, 0.042, 1)
    br.inputs['Color2'].default_value = (0.105, 0.048, 0.034, 1)
    br.inputs['Mortar'].default_value = (0.135, 0.130, 0.118, 1)
    br.inputs['Mortar Size'].default_value = 0.022
    nt.links.new(tc.outputs['Object'], br.inputs['Vector'])
    n = noise_node(nt, (-600, -300), 4.0, 7.0, 0.6, links_from=tc.outputs['Object'])
    mx = mix_node(nt, (-330, 0))
    mix_set(mx, b=(0.060, 0.085, 0.030, 1), fac=0.35)
    nt.links.new(br.outputs['Color'], mx.inputs[6])
    nt.links.new(n.outputs['Fac'], mx.inputs[0])
    nt.links.new(mix_out(mx), b.inputs['Base Color'])
    sset(b, 'Roughness', 0.95)
    bump = nt.nodes.new('ShaderNodeBump'); bump.location = (120, -300)
    sset(bump, 'Strength', 0.5)
    nt.links.new(br.outputs['Fac'], bump.inputs['Height'])
    nt.links.new(bump.outputs['Normal'], b.inputs['Normal'])
    return m

M_RUST = mat_rust()
M_WALL = mat_wall()
M_WOOD = mat_wood()
M_LEAF = mat_leaf()
M_GRASS = mat_grass()
M_GROUND = mat_ground()
M_STONE = mat_stone()
M_BRICK = mat_brick()

bms = {k: bmesh.new() for k in ('wall', 'roof', 'wood', 'branch', 'leaf', 'stone', 'brick')}

# ----------------------------------------------------------------------------- ground
def gh(x, y):
    v = 0.34 * mnoise.noise(Vector((x * 0.045, y * 0.045, 0.5)))
    v += 0.11 * mnoise.noise(Vector((x * 0.16, y * 0.16, 1.7)))
    pathf = max(0.0, 1.0 - abs(x) / 3.2)
    return v * (1.0 - 0.75 * pathf) - 0.07 * pathf

GN, GS = 120, 38.0
gbm = bmesh.new()
grid = []
for i in range(GN + 1):
    row = []
    x = -GS + 2 * GS * i / GN
    for j in range(GN + 1):
        y = -GS + 2 * GS * j / GN
        row.append(gbm.verts.new((x, y, gh(x, y))))
    grid.append(row)
for i in range(GN):
    for j in range(GN):
        try:
            gbm.faces.new([grid[i][j], grid[i + 1][j], grid[i + 1][j + 1], grid[i][j + 1]])
        except Exception:
            pass
obj_from_bm("Ground", gbm, M_GROUND, smooth=True)

# ----------------------------------------------------------------------------- house
def build_wall(Mw, length, height, thickness, openings, rnd, outward):
    ops = sorted(openings, key=lambda o: o[0] - o[1] * 0.5)
    cursor = -length * 0.5
    for (cx, ow, z0, z1) in ops:
        x0, x1 = cx - ow * 0.5, cx + ow * 0.5
        x0 = max(x0, -length * 0.5); x1 = min(x1, length * 0.5)
        if x0 > cursor + 1e-4:
            wdt = x0 - cursor
            add_box(bms['wall'], Mw @ Matrix.Translation(((cursor + x0) * 0.5, 0, height * 0.5)),
                    (wdt, thickness, height))
        if z0 > 1e-4:
            add_box(bms['wall'], Mw @ Matrix.Translation(((x0 + x1) * 0.5, 0, z0 * 0.5)),
                    (x1 - x0, thickness, z0))
        if z1 < height - 1e-4:
            add_box(bms['wall'], Mw @ Matrix.Translation(((x0 + x1) * 0.5, 0, (z1 + height) * 0.5)),
                    (x1 - x0, thickness, height - z1))
        cursor = x1
    if cursor < length * 0.5 - 1e-4:
        wdt = length * 0.5 - cursor
        add_box(bms['wall'], Mw @ Matrix.Translation(((cursor + length * 0.5) * 0.5, 0, height * 0.5)),
                (wdt, thickness, height))
    # frames, shutters, foliage bursting out
    for (cx, ow, z0, z1) in ops:
        if z1 - z0 < 0.5 or ow < 0.4:
            continue
        if z1 > height - 0.05:
            continue
        add_box(bms['wood'], Mw @ Matrix.Translation((cx, 0, z0 - 0.04)),
                (ow + 0.20, thickness * 1.5, 0.10))
        add_box(bms['wood'], Mw @ Matrix.Translation((cx, 0, z1 + 0.05)),
                (ow + 0.24, thickness * 1.6, 0.12))
        for s in (-1, 1):
            add_box(bms['wood'], Mw @ Matrix.Translation((cx + s * (ow * 0.5 + 0.05), 0, (z0 + z1) * 0.5)),
                    (0.10, thickness * 1.4, z1 - z0))
        if rnd.random() < 0.65:
            ang = rnd.uniform(0.15, 0.6) * rnd.choice((-1, 1))
            sM = (Mw @ Matrix.Translation((cx + rnd.choice((-1, 1)) * (ow * 0.5 + 0.22), -0.14, (z0 + z1) * 0.5))
                  @ Matrix.Rotation(ang, 4, 'Y'))
            add_box(bms['wood'], sM, (ow * 0.55, 0.05, (z1 - z0) * 0.9))
        c = Mw @ Vector((cx, 0.0, (z0 + z1) * 0.5))
        leaf_cluster(bms['leaf'], Vector(c) + outward * 0.16,
                     min(ow, z1 - z0) * 0.62, 170, 0.115, rnd, flatten=0.9)

def add_gable(bm, M, x, thickness, d, h, rise):
    pts = [(-d * 0.5, h), (d * 0.5, h), (0.0, h + rise)]
    layers = []
    for s in (-0.5, 0.5):
        layers.append([bm.verts.new(M @ Vector((x + s * thickness, y, z))) for (y, z) in pts])
    try:
        bm.faces.new(layers[0])
        bm.faces.new(layers[1][::-1])
    except Exception:
        pass
    for i in range(3):
        j = (i + 1) % 3
        try:
            bm.faces.new([layers[0][i], layers[0][j], layers[1][j], layers[1][i]])
        except Exception:
            pass

def add_vine(M, base_local, up, side, length, rnd, leafsize=0.13):
    pts = []
    n = 16
    ph = rnd.uniform(0, 3.0)
    fr = rnd.uniform(3.5, 7.0)
    for i in range(n + 1):
        t = i / n
        off = math.sin(t * fr + ph) * 0.42 * (0.25 + t)
        pts.append(base_local + up * (length * t) + side * off)
    add_tube(bms['branch'], pts, 0.028, sides=5)
    for i, p in enumerate(pts):
        if i % 2 == 0:
            leaf_cluster(bms['leaf'], p, 0.22, 22, leafsize, rnd, flatten=1.0)

roof_spots = []
blockers = []

def build_house(loc, rotz, w, d, h, rise=1.55, ov=0.5, seed=1, ruin=False, chimney=True):
    rnd = random.Random(seed)
    z = gh(loc[0], loc[1])
    M = (Matrix.Translation(Vector((loc[0], loc[1], z)))
         @ Euler((rnd.uniform(-0.012, 0.012), rnd.uniform(-0.018, 0.018), rotz), 'XYZ').to_matrix().to_4x4())
    R3 = M.to_3x3()
    t = 0.22
    add_box(bms['stone'], M @ Matrix.Translation((0, 0, -0.32)), (w + 0.55, d + 0.55, 0.75))

    # front
    fw = M @ Matrix.Translation((0, -d * 0.5, 0))
    ops_f = [(-w * 0.20, 1.05, 0.0, 2.15), (w * 0.27, 1.05, 1.15, 2.25)]
    if w > 6.0:
        ops_f.append((-w * 0.42, 0.9, 1.2, 2.2))
    build_wall(fw, w, h, t, ops_f, rnd, (R3 @ Vector((0, -1, 0))).normalized())
    # back
    bw = M @ Matrix.Translation((0, d * 0.5, 0))
    ops_b = [(rnd.uniform(-0.6, 0.6), 1.15, 1.15, 2.3)]
    if ruin:
        ops_b.append((w * 0.30, 2.3, 1.7, h + 0.1))
    build_wall(bw, w, h, t, ops_b, rnd, (R3 @ Vector((0, 1, 0))).normalized())
    # sides
    for sgn in (-1, 1):
        sw = M @ Matrix.Translation((sgn * (w * 0.5), 0, 0)) @ Matrix.Rotation(math.pi * 0.5, 4, 'Z')
        ops_s = [(rnd.uniform(-0.7, 0.7), 1.0, 1.25, 2.25)]
        build_wall(sw, d, h, t, ops_s, rnd, (R3 @ Vector((sgn, 0, 0))).normalized())
        add_gable(bms['wall'], M, sgn * (w * 0.5), t, d, h, rise)

    # roof
    p = math.atan2(rise, d * 0.5 + ov)
    ly = math.sqrt(rise ** 2 + (d * 0.5 + ov) ** 2)
    total = w + 2 * ov
    nstrip = max(4, int(round(total / 0.95)))
    sw_ = total / nstrip
    for side in (-1, 1):
        for k in range(nstrip):
            cx = -total * 0.5 + sw_ * (k + 0.5)
            r = rnd.random()
            if r < 0.10 and not (k in (0, nstrip - 1)):
                continue  # missing panel -> hole
            y1 = ly
            y0 = 0.0
            if r < 0.20:
                y1 = ly * rnd.uniform(0.45, 0.8)   # torn / partly gone
            base = M @ Matrix.Translation((cx, 0, h + rise))
            if side > 0:
                Mr = base @ Matrix.Rotation(-p, 4, 'X')
            else:
                Mr = base @ Matrix.Rotation(math.pi, 4, 'Z') @ Matrix.Rotation(-p, 4, 'X')
            add_corrugated(bms['roof'], Mr, sw_ * 0.99, y1, amp=0.035, period=0.28,
                           nv=5, sag=rnd.uniform(0.0, 0.08), y0=y0)
            # moss/grass on the roof plane
            for _ in range(3):
                v = rnd.uniform(0.15, min(0.95, y1 / ly))
                wp = Mr @ Vector((rnd.uniform(-sw_ * 0.4, sw_ * 0.4), v * ly, 0.05))
                if rnd.random() < 0.6:
                    leaf_cluster(bms['leaf'], wp, 0.28, 40, 0.09, rnd, flatten=0.35)
                else:
                    roof_spots.append((Mr, rnd.uniform(-sw_ * 0.4, sw_ * 0.4), v * ly, side))
    # ridge cap
    add_box(bms['roof'], M @ Matrix.Translation((0, 0, h + rise + 0.05)), (total * 0.98, 0.34, 0.10))
    # purlins & ceiling joists visible through holes
    for side in (-1, 1):
        for v in (0.30, 0.68):
            yy = side * v * (d * 0.5 + ov)
            zz = h + rise - v * rise - 0.10
            add_box(bms['wood'], M @ Matrix.Translation((0, yy, zz)), (w + 0.3, 0.10, 0.14))
    nj = max(3, int(w / 1.4))
    for k in range(nj):
        xx = -w * 0.5 + w * (k + 0.5) / nj
        add_box(bms['wood'], M @ Matrix.Translation((xx, 0, h - 0.10)), (0.10, d, 0.16))
    # chimney
    if chimney:
        cxx = w * rnd.choice((-0.28, 0.30))
        add_box(bms['brick'], M @ Matrix.Translation((cxx, 0.0, (h + rise + 0.75) * 0.5)),
                (0.75, 0.75, h + rise + 0.75))
        leaf_cluster(bms['leaf'], M @ Vector((cxx, 0, h + rise + 0.7)), 0.4, 45, 0.1, rnd, flatten=0.5)
    # rubble at ruined wall
    if ruin:
        for _ in range(26):
            pl = Vector((rnd.uniform(-w * 0.45, w * 0.45), d * 0.5 + rnd.uniform(0.3, 1.6),
                         rnd.uniform(0.0, 0.35)))
            s = rnd.uniform(0.14, 0.42)
            add_box(bms['stone'], M @ Matrix.Translation(pl)
                    @ Euler((rnd.uniform(0, 3), rnd.uniform(0, 3), rnd.uniform(0, 3))).to_matrix().to_4x4(),
                    (s, s * rnd.uniform(0.6, 1.4), s * rnd.uniform(0.4, 0.9)))
    # ivy up the walls
    up = Vector((0, 0, 1))
    for _ in range(3):
        wchoice = rnd.choice(('f', 'b', 'l', 'r'))
        if wchoice == 'f':
            base = Vector((rnd.uniform(-w * 0.45, w * 0.45), -d * 0.5 - t * 0.6 - 0.06, 0.02))
            side = Vector((1, 0, 0))
        elif wchoice == 'b':
            base = Vector((rnd.uniform(-w * 0.45, w * 0.45), d * 0.5 + t * 0.6 + 0.06, 0.02))
            side = Vector((1, 0, 0))
        elif wchoice == 'l':
            base = Vector((-w * 0.5 - t * 0.6 - 0.06, rnd.uniform(-d * 0.4, d * 0.4), 0.02))
            side = Vector((0, 1, 0))
        else:
            base = Vector((w * 0.5 + t * 0.6 + 0.06, rnd.uniform(-d * 0.4, d * 0.4), 0.02))
            side = Vector((0, 1, 0))
        gm = M
        pts_base = gm @ base
        u = (R3 @ up).normalized()
        s2 = (R3 @ side).normalized()
        add_vine(gm, Vector(pts_base), u, s2, h * rnd.uniform(0.75, 1.25), rnd)
    # bushes hugging the base
    for _ in range(10):
        e = rnd.choice(('f', 'b', 'l', 'r'))
        if e == 'f':
            bl = Vector((rnd.uniform(-w * 0.5, w * 0.5), -d * 0.5 - rnd.uniform(0.15, 0.7), 0.05))
        elif e == 'b':
            bl = Vector((rnd.uniform(-w * 0.5, w * 0.5), d * 0.5 + rnd.uniform(0.15, 0.7), 0.05))
        elif e == 'l':
            bl = Vector((-w * 0.5 - rnd.uniform(0.15, 0.7), rnd.uniform(-d * 0.5, d * 0.5), 0.05))
        else:
            bl = Vector((w * 0.5 + rnd.uniform(0.15, 0.7), rnd.uniform(-d * 0.5, d * 0.5), 0.05))
        leaf_cluster(bms['leaf'], M @ bl, rnd.uniform(0.35, 0.75), 140, 0.11, rnd, flatten=0.8)
    blockers.append((loc[0], loc[1], max(w, d) * 0.62 + 0.9))
    return M

HOUSES = [
    dict(loc=(-6.9, 0.5), rotz=0.06, w=6.2, d=5.2, h=3.2, seed=3, ruin=False),
    dict(loc=(7.4, 5.0), rotz=math.pi - 0.09, w=7.0, d=5.6, h=3.5, seed=8, ruin=True),
    dict(loc=(-7.9, 11.8), rotz=0.14, w=5.8, d=5.0, h=3.0, seed=15, ruin=True),
    dict(loc=(8.1, -4.2), rotz=math.pi + 0.05, w=5.2, d=4.6, h=2.9, seed=21, ruin=False),
    dict(loc=(-9.0, 21.5), rotz=-0.10, w=5.0, d=4.6, h=2.9, seed=29, ruin=False),
    dict(loc=(7.9, 16.5), rotz=math.pi + 0.35, w=3.6, d=3.2, h=2.3, seed=34, ruin=True, chimney=False),
]
house_mats = []
for hh in HOUSES:
    house_mats.append(build_house(hh['loc'], hh['rotz'], hh['w'], hh['d'], hh['h'],
                                  seed=hh['seed'], ruin=hh['ruin'],
                                  chimney=hh.get('chimney', True)))

# ----------------------------------------------------------------------------- trees
def build_tree(x, y, height, seed, dead=False):
    rnd = random.Random(seed)
    base = Vector((x, y, gh(x, y) - 0.1))
    n = 9
    pts = []
    radii = []
    r0 = height * 0.055
    ox = oy = 0.0
    for i in range(n + 1):
        t = i / n
        pts.append(base + Vector((ox, oy, height * 0.55 * t)))
        radii.append(r0 * (1 - 0.62 * t))
        ox += rnd.uniform(-0.09, 0.09)
        oy += rnd.uniform(-0.09, 0.09)
    add_tube(bms['branch'], pts, 0, sides=8, radii=radii)
    tips = []
    for k in range(6):
        i0 = int(n * rnd.uniform(0.4, 0.98))
        start = pts[i0]
        ang = rnd.uniform(0, TAU)
        L = height * rnd.uniform(0.32, 0.55)
        d = Vector((math.cos(ang), math.sin(ang), rnd.uniform(0.7, 1.6))).normalized()
        bp, br = [], []
        m = 5
        for j in range(m + 1):
            tt = j / m
            p = (start + d * (L * tt) + Vector((0, 0, -0.22 * L * tt * tt))
                 + Vector((rnd.uniform(-0.12, 0.12), rnd.uniform(-0.12, 0.12), 0)))
            bp.append(p)
            br.append(radii[i0] * 0.7 * (1 - 0.8 * tt) + 0.012)
        add_tube(bms['branch'], bp, 0, sides=6, radii=br)
        tips.append(bp[-1])
    if not dead:
        for tp in tips:
            leaf_cluster(bms['leaf'], tp, height * 0.17, 300, height * 0.030, rnd, flatten=0.85)
        leaf_cluster(bms['leaf'], pts[-1] + Vector((0, 0, height * 0.08)),
                     height * 0.16, 220, height * 0.030, rnd, flatten=0.8)
    # ivy on trunk
    add_vine(Matrix(), base + Vector((r0, 0, 0.05)), Vector((0, 0, 1)), Vector((1, 0, 0)),
             height * 0.4, rnd, leafsize=0.1)
    blockers.append((x, y, 1.5))

build_tree(-14.5, 15.0, 9.5, 101)
build_tree(12.5, 19.0, 8.5, 202)
build_tree(-15.0, -3.5, 10.5, 303)
build_tree(14.0, -9.0, 8.0, 404, dead=True)
build_tree(-11.0, 27.5, 8.0, 505)
build_tree(11.5, 28.0, 9.0, 606)
build_tree(-4.5, 33.0, 9.5, 707)
build_tree(4.8, 34.0, 8.5, 808, dead=True)

# ----------------------------------------------------------------------------- fences, debris
rnd = random.Random(77)
for side in (-1, 1):
    y = -12.0
    while y < 30.0:
        x = side * (2.5 + rnd.uniform(-0.2, 0.2))
        if random.random() < 0.25:
            y += rnd.uniform(1.2, 2.6)
            continue
        z = gh(x, y)
        Mp = (Matrix.Translation(Vector((x, y, z + 0.45)))
              @ Euler((rnd.uniform(-0.25, 0.25), rnd.uniform(-0.25, 0.25), rnd.uniform(0, 1))).to_matrix().to_4x4())
        add_box(bms['wood'], Mp, (0.12, 0.12, rnd.uniform(0.7, 1.15)))
        if rnd.random() < 0.5:
            leaf_cluster(bms['leaf'], Vector((x, y, z + 0.25)), 0.35, 60, 0.1, rnd, flatten=0.7)
        y += rnd.uniform(1.4, 2.4)

# fallen corrugated sheets
for (fx, fy, fr) in [(-3.6, -1.5, 0.5), (4.2, 2.6, -0.7), (-4.4, 9.5, 1.4), (5.1, 12.5, 0.3),
                     (-3.0, 18.0, -1.1)]:
    z = gh(fx, fy)
    Ms = (Matrix.Translation(Vector((fx, fy, z + 0.08)))
          @ Euler((rnd.uniform(-0.12, 0.12), rnd.uniform(-0.12, 0.12), fr)).to_matrix().to_4x4())
    add_corrugated(bms['roof'], Ms, rnd.uniform(0.9, 1.5), rnd.uniform(1.6, 2.6),
                   amp=0.04, period=0.28, nv=4, sag=0.06)
    leaf_cluster(bms['leaf'], Vector((fx, fy, z + 0.1)), 0.7, 90, 0.1, rnd, flatten=0.35)

# barrel / crate debris
for (bx, by) in [(-4.0, 3.5), (5.6, -1.4), (-5.2, 15.0), (6.4, 9.4)]:
    z = gh(bx, by)
    add_box(bms['wood'], Matrix.Translation(Vector((bx, by, z + 0.3)))
            @ Euler((rnd.uniform(-0.2, 0.2), rnd.uniform(-0.2, 0.2), rnd.uniform(0, 2))).to_matrix().to_4x4(),
            (0.7, 0.7, 0.65))
    leaf_cluster(bms['leaf'], Vector((bx, by, z + 0.55)), 0.45, 70, 0.1, rnd, flatten=0.5)

# ----------------------------------------------------------------------------- realise merged meshes
obj_from_bm("Walls", bms['wall'], M_WALL)
obj_from_bm("RustRoofs", bms['roof'], M_RUST, smooth=True)
obj_from_bm("Timber", bms['wood'], M_WOOD)
obj_from_bm("Branches", bms['branch'], M_WOOD, smooth=True)
obj_from_bm("Foliage", bms['leaf'], M_LEAF)
obj_from_bm("Rubble", bms['stone'], M_STONE)
obj_from_bm("Brickwork", bms['brick'], M_BRICK)

# ----------------------------------------------------------------------------- scattered instances
def make_tuft(name, nblades, hmin, hmax):
    bm = bmesh.new()
    r = random.Random(hash(name) & 0xffff)
    for _ in range(nblades):
        a = r.uniform(0, TAU)
        rr = r.uniform(0, 0.11)
        base = Vector((math.cos(a) * rr, math.sin(a) * rr, 0))
        add_blade(bm, base, r.uniform(hmin, hmax), r.uniform(0, TAU), r.uniform(0.05, 0.25), 0.022)
    me = bpy.data.meshes.new(name)
    bm.to_mesh(me); bm.free()
    me.materials.append(M_GRASS)
    return me

TUFT = make_tuft("GrassTuft", 16, 0.20, 0.48)
TUFT_S = make_tuft("GrassTuftSmall", 9, 0.10, 0.24)

def make_bush(name, radius, count, size, flat):
    bm = bmesh.new()
    r = random.Random(hash(name) & 0xffff)
    leaf_cluster(bm, Vector((0, 0, radius * 0.55 * flat)), radius, count, size, r, flatten=flat)
    me = bpy.data.meshes.new(name)
    bm.to_mesh(me); bm.free()
    me.materials.append(M_LEAF)
    return me

BUSH = make_bush("Bush", 0.7, 300, 0.12, 0.8)
BUSH_S = make_bush("BushSmall", 0.35, 120, 0.10, 0.9)

def make_rock(name):
    bm = bmesh.new()
    bmesh.ops.create_icosphere(bm, subdivisions=2, radius=1.0)
    r = random.Random(hash(name) & 0xffff)
    for v in bm.verts:
        v.co += Vector((r.uniform(-0.22, 0.22), r.uniform(-0.22, 0.22), r.uniform(-0.18, 0.18)))
        v.co.z *= 0.6
    me = bpy.data.meshes.new(name)
    bm.to_mesh(me); bm.free()
    me.materials.append(M_STONE)
    return me

ROCK = make_rock("Rock")

def instance(me, name, loc, rotz, scale, matrix=None):
    ob = bpy.data.objects.new(name, me)
    if matrix is not None:
        ob.matrix_world = matrix
    else:
        ob.location = loc
        ob.rotation_euler = (0, 0, rotz)
        ob.scale = (scale, scale, scale)
    return link(ob)

def free_spot(x, y, pad=0.0):
    for (hx, hy, r) in blockers:
        if (x - hx) ** 2 + (y - hy) ** 2 < (r + pad) ** 2:
            return False
    return True

placed = 0
tries = 0
srnd = random.Random(4242)
while placed < 1100 and tries < 12000:
    tries += 1
    x = srnd.uniform(-32, 32)
    y = srnd.uniform(-20, 35)
    if not free_spot(x, y):
        continue
    on_path = abs(x) < 2.1
    if on_path and srnd.random() < 0.9:
        continue
    z = gh(x, y) - 0.03
    me = TUFT_S if on_path or srnd.random() < 0.3 else TUFT
    sc = srnd.uniform(0.6, 1.6) * (0.6 if on_path else 1.0)
    instance(me, "Grass", (x, y, z), srnd.uniform(0, TAU), sc)
    placed += 1

for _ in range(90):
    for _t in range(20):
        x = srnd.uniform(-30, 30); y = srnd.uniform(-18, 33)
        if free_spot(x, y, 0.3) and abs(x) > 2.4:
            break
    z = gh(x, y)
    me = BUSH if srnd.random() < 0.55 else BUSH_S
    instance(me, "Bush", (x, y, z), srnd.uniform(0, TAU), srnd.uniform(0.7, 1.7))

for _ in range(70):
    x = srnd.uniform(-30, 30); y = srnd.uniform(-18, 33)
    z = gh(x, y)
    s = srnd.uniform(0.12, 0.4)
    ob = instance(ROCK, "Rock", (x, y, z + s * 0.25), srnd.uniform(0, TAU), 1.0)
    ob.scale = (s * srnd.uniform(0.8, 1.5), s * srnd.uniform(0.8, 1.5), s)

# grass on rooftops (aligned to the slope)
rr = random.Random(999)
for (Mr, lx, lyv, side) in roof_spots:
    Mfull = Mr @ Matrix.Translation((lx, lyv, 0.03)) @ Matrix.Diagonal(
        Vector((rr.uniform(0.6, 1.2),) * 3 + (1.0,)))
    instance(TUFT_S, "RoofGrass", (0, 0, 0), 0, 1.0, matrix=Mfull)

# ----------------------------------------------------------------------------- lighting / world
world = bpy.data.worlds.new("VillageWorld")
scene.world = world
world.use_nodes = True
wnt = world.node_tree
for n in list(wnt.nodes):
    wnt.nodes.remove(n)
wout = wnt.nodes.new('ShaderNodeOutputWorld'); wout.location = (400, 0)
bg = wnt.nodes.new('ShaderNodeBackground'); bg.location = (180, 0)
sky = wnt.nodes.new('ShaderNodeTexSky'); sky.location = (-100, 0)
try:
    sky.sky_type = 'NISHITA'
    sky.sun_elevation = math.radians(21.0)
    sky.sun_rotation = math.radians(232.0)
    sky.altitude = 200.0
    sky.air_density = 1.5
    sky.dust_density = 3.2
    sky.sun_intensity = 0.5
    sky.sun_disc = False
except Exception:
    pass
bg.inputs['Strength'].default_value = 0.85
wnt.links.new(sky.outputs[0], bg.inputs['Color'])
wnt.links.new(bg.outputs[0], wout.inputs['Surface'])

sun_data = bpy.data.lights.new("Sun", type='SUN')
sun_data.energy = 3.4
sun_data.color = (1.0, 0.83, 0.60)
sun_data.angle = math.radians(1.8)
sun = link(bpy.data.objects.new("Sun", sun_data))
sun.rotation_euler = Euler((math.radians(69.0), 0.0, math.radians(232.0)), 'XYZ')

fill_data = bpy.data.lights.new("SkyFill", type='AREA')
fill_data.energy = 900.0
fill_data.size = 22.0
fill_data.color = (0.62, 0.74, 0.95)
fill = link(bpy.data.objects.new("SkyFill", fill_data))
fill.location = (-16.0, -14.0, 16.0)
fill.rotation_euler = (math.radians(52), 0, math.radians(-42))

# ----------------------------------------------------------------------------- camera
cam_data = bpy.data.cameras.new("Camera")
cam_data.lens = 38.0
cam_data.clip_end = 400.0
cam = link(bpy.data.objects.new("Camera", cam_data))
cam.location = Vector((2.0, -14.8, 3.7))
target = Vector((-0.8, 5.0, 1.9))
d = target - cam.location
cam.rotation_euler = d.to_track_quat('-Z', 'Y').to_euler()
scene.camera = cam

# ----------------------------------------------------------------------------- render settings
scene.render.engine = 'CYCLES'
scene.cycles.samples = 64
scene.cycles.use_denoising = True
scene.cycles.max_bounces = 8
scene.cycles.transmission_bounces = 4
scene.cycles.transparent_max_bounces = 8
try:
    scene.cycles.device = 'GPU'
    prefs = bpy.context.preferences.addons.get('cycles')
    if prefs:
        cp = prefs.preferences
        for t in ('OPTIX', 'CUDA', 'HIP', 'METAL', 'ONEAPI'):
            try:
                cp.compute_device_type = t
                break
            except Exception:
                continue
        try:
            cp.get_devices()
            for dev in cp.devices:
                dev.use = True
        except Exception:
            pass
except Exception:
    pass
scene.render.resolution_x = 1536
scene.render.resolution_y = 1536
scene.render.film_transparent = False
try:
    scene.view_settings.view_transform = 'AgX'
except Exception:
    pass
try:
    scene.view_settings.look = 'AgX - Medium Contrast'
except Exception:
    pass
scene.view_settings.exposure = 0.15

try:
    if bpy.data.filepath:
        bpy.ops.wm.save_mainfile()
except Exception:
    pass
```