Program za generisanje i sečenje 3V geodetske kupole

make_dome_clip_yrotate.py je Python skript namenjen za generisanje 3V geodetske kupole, njenu rotaciju i sečenje, i na kraju izvoz u više CAD i 3D formata. Posebno je prilagođen za slučaj gde se dužine letvica A, B i C unapred dobijaju iz drugog programa, geosf.py koji sračunava dužinu, broj letvica i broj potrebnih konektora za izgradnju Geodetske kupole pa se kupola konstruiše na osnovu tih realnih vrednosti.

Osnovne funkcije programa

1. Generisanje geodetske strukture
   • Polazna figura je ikosaedar (20 jednakostraničnih trouglova).
   • Svaka stranica se deli na manje delove (frekvencija freq=3, tj. 3V), a zatim se svi dobijeni čvorovi projektuju na sferu poluprečnika R=500 cm.
   • Na taj način se dobija mreža geodetske kupole sa različitim tipovima ivica.
2. Rotacija kupole
   • Ceo sistem čvorova se rotira oko Y ose za 58.2825°.
   • Ova rotacija nije proizvoljna – ona omogućava da kupola nakon presecanja dobije željeni oblik i proporcije.
3. Odsecanje kupole
   • Nakon rotacije, primenjuje se sečenje ravni z = z_clip.
   • Visina presecanja zavisi od parametra CLIP_MODE:
     • "3/8" – odsecanje na 1/4 radijusa iznad centra.
     • "5/8" – odsecanje na 1/4 radijusa ispod centra.
     • "manual" – korisnik ručno unosi visinu presecanja.
   • Time se dobija samo gornji deo kupole, kao što je slučaj kod kupola koje se u praksi oslanjaju na horizontalnu osnovu.
4. Dužine letvica i klasifikacija
   • Svaka ivica se meri i poredi sa unapred datim vrednostima:
     • A = 174.31 cm
     • B = 201.78 cm
     • C = 206.21 cm
   • Program svaku ivicu svrstava u tip A, B ili C, u skladu sa najbližom dužinom.
   • Na kraju se ispisuje broj elemenata po tipovima, što omogućava proveru da li se dobijaju očekivane količine (30 A, 40 B, 50 C).
5. Modelovanje strutova (letvica)
   • Svaka ivica se modeluje kao cilindar određenog radijusa (2.5 cm).
   • Cilindri se generišu pomoću biblioteke trimesh, i seku se na ravni z_clip kako bi samo deo iznad preseka ostao u modelu.
6. Izvoz u formate
   Program rezultate snima u tri formata:
   • PLY (dome_clipped_y_colored.ply) – sa bojama po tipu letvica (A=crvena, B=zelena, C=plava).
   • STL (dome_clipped_y_combined.stl) – kombinovani 3D model spreman za štampu ili dalju obradu.
   • DXF (dome_clipped_y_centerlines.dxf) – nacrt sa centar-linijama svake letvice, grupisanim u slojeve po tipovima (A/B/C) i obojenim prema standardnim CAD bojama.

Tehnički detalji

• Ulazni parametri:
  • poluprečnik kupole (R),
  • frekvencija subdivizije (freq),
  • dužine letvica (iz geosf.py programa),
  • način presecanja (CLIP_MODE).
• Izlazni fajlovi:
  • PLY sa bojama,
  • STL za 3D modeliranje i štampu,
  • DXF za tehničku dokumentaciju i proizvodnju.
• Jedinice: svi proračuni i fajlovi su u centimetrima.

Primena

Ovaj program je koristan u:

• projektovanju geodetskih kupola za arhitektonske i inženjerske svrhe,
• pripremi dokumentacije za izradu konstrukcije od letvica,
• analizi količina materijala i tipova letvica,
• izvozu podataka u 3D i CAD programe (npr. AutoCAD, Blender, Rhino).

📌 Ukratko: make_dome_clip_yrotate.py je specijalizovan alat za projektovanje 3V geodetskih kupola, koji omogućava realističan proračun dužina letvica, njihovu klasifikaciju, sečenje kupole i izvoz u standardne tehničke formate.

# The MIT License (MIT)
# Copyright (c) 2025 Aleksandar Maričić
#
# Ovim se omogućava bilo kome da koristi, kopira, menja, spaja, objavljuje,
# distribuira, daje podlicencu i/ili prodaje kopije ovog softverskog programa,
# uz uslov da u svim kopijama ili značajnim delovima softverskog programa bude
# uključena sledeća obavest:
#
# Copyright (c) 2025 Aleksandar Maričić
#
# Ovaj softverski program je pružen "takav kakav jeste", bez bilo kakvih garancija,
# izričitih ili impliciranih, uključujući, ali ne ograničavajući se na, garancije o
# prikladnosti za prodaju ili pogodnosti za određenu svrhu. U svakom slučaju, autori
# ili nosioci prava nisu odgovorni za bilo kakvu štetu ili druge obaveze koje mogu nastati
# usled upotrebe ovog softverskog programa.





#!/usr/bin/env python3
"""
make_dome_clip_yrotate.py

Generiše 3V geodetsku kupolu (freq=3), rotira oko Y ose za 58.2825°,
odseca sve ispod ravni z=Z_CLIP_HEIGHT (posle rotacije).
Snima PLY (boje po tipu A/B/C), STL i DXF (centerlines).

Jedinice: cm
"""
import numpy as np
import math
import trimesh
import ezdxf
from collections import Counter

# ---------- PODESAVANJA ----------
R = 500.0                # poluprečnik u cm
freq = 3                 # 3V
strut_radius = 2.5       # cm - poluprečnik cilindra
cylinder_sections = 20
ROT_Y_DEG = 58.2825      # stepeni rotacije oko Y ose
GIVEN = {'A': 174.31, 'B': 201.78, 'C': 206.21}
EXPECTED_COUNTS = {'A': 30, 'B': 40, 'C': 50}
OUT_PLY = "dome_clipped_y_colored.ply"
OUT_STL = "dome_clipped_y_combined.stl"
OUT_DXF = "dome_clipped_y_centerlines.dxf"
Z_TOL = 1e-8             # tolerancija pri sečenju

# --- kontrola presecanja ---
CLIP_MODE = "5/8"        # "3/8", "5/8", "manual"
Z_CLIP_HEIGHT_MANUAL = 0.1
# ---------------------------------

def compute_z_clip(R, mode, manual_value):
    if mode == "3/8":
        return R/4.0
    elif mode == "5/8":
        return -R/4.0
    elif mode == "manual":
        return manual_value
    else:
        raise ValueError(f"Nepoznat CLIP_MODE: {mode}")

# ----- generisanje icosahedra -----
def make_icosahedron_unit():
    phi = (1.0 + 5.0**0.5) / 2.0
    V = np.array([
        (-1,  phi,  0),( 1,  phi,  0),(-1, -phi,  0),( 1, -phi,  0),
        ( 0, -1,  phi),( 0,  1,  phi),( 0, -1, -phi),( 0,  1, -phi),
        ( phi,  0, -1),( phi,  0,  1),(-phi,  0, -1),(-phi,  0,  1)
    ], dtype=float)
    V = V / np.linalg.norm(V, axis=1)[:,None]
    F = np.array([
    (0,11,5),(0,5,1),(0,1,7),(0,7,10),(0,10,11),
    (1,5,9),(5,11,4),(11,10,2),(10,7,6),(7,1,8),
    (3,9,4),(3,4,2),(3,2,6),(3,6,8),(3,8,9),
    (4,9,5),(2,4,11),(6,2,10),(8,6,7),(9,8,1)
    ])
    return V, F

def subdivide_and_project(V, F, freq, R):
    pts = []
    faces = []
    pt_map = {}
    def key_of(p):
        return (round(float(p[0]),8), round(float(p[1]),8), round(float(p[2]),8))
    for tri in F:
        v0, v1, v2 = V[list(tri)]
        local_idx = {}
        for i in range(freq+1):
            for j in range(freq+1-i):
                k = freq - i - j
                p = (i * v0 + j * v1 + k * v2) / freq
                p = p / np.linalg.norm(p) * R
                key = key_of(p)
                if key not in pt_map:
                    pt_map[key] = len(pts)
                    pts.append(p)
                local_idx[(i,j)] = pt_map[key]
        for i in range(freq):
            for j in range(freq - i):
                a = local_idx[(i,j)]
                b = local_idx[(i+1,j)]
                c = local_idx[(i,j+1)]
                faces.append((a,b,c))
                if j + i + 1 < freq:
                    d = local_idx[(i+1,j+1)]
                    faces.append((b,d,c))
    return np.array(pts), np.array(faces)

def build_edges(faces):
    edges = {}
    for (a,b,c) in faces:
        for u,v in ((a,b),(b,c),(c,a)):
            e = tuple(sorted((int(u),int(v))))
            edges[e] = edges.get(e, 0) + 1
    return list(edges.keys())

def edge_length(pts, e):
    u,v = e
    return np.linalg.norm(pts[u]-pts[v])

def nearest_label(length, given):
    best = None
    bestd = 1e9
    for label,val in given.items():
        d = abs(length - val)
        if d < bestd:
            bestd = d; best = label
    return best

def rotation_matrix_y(deg):
    th = math.radians(deg)
    c = math.cos(th); s = math.sin(th)
    Rm = np.array([[c,0,s],[0,1,0],[-s,0,c]], dtype=float)
    return Rm

def cylinder_between_clipped(p0, p1, radius, sections, z_clip):
    height = np.linalg.norm(p1 - p0)
    if height <= 1e-9:
        return None
    cyl = trimesh.creation.cylinder(radius=radius, height=height, sections=sections)
    mid = (p0 + p1) / 2.0
    direction = (p1 - p0) / height
    z = np.array([0.0, 0.0, 1.0])
    v = np.cross(z, direction)
    c = np.dot(z, direction)
    if np.linalg.norm(v) < 1e-8 and c > 0.999999:
        Rm = np.eye(3)
    elif np.linalg.norm(v) < 1e-8 and c < -0.999999:
        Rm = trimesh.transformations.rotation_matrix(math.pi, [1,0,0])[:3,:3]
    else:
        s = np.linalg.norm(v)
        kmat = np.array([[0,-v[2],v[1]],[v[2],0,-v[0]],[-v[1],v[0],0]])
        Rm = np.eye(3) + kmat + (kmat @ kmat) * ((1 - c) / (s**2))
    transform = np.eye(4)
    transform[:3,:3] = Rm
    transform[:3,3] = mid
    cyl.apply_transform(transform)

    centroids = cyl.triangles_center
    mask = centroids[:,2] >= z_clip - Z_TOL
    if not mask.any():
        return None
    faces_idx = np.nonzero(mask)[0]
    faces = cyl.faces[faces_idx]
    unique_vid, inverse = np.unique(faces.flatten(), return_inverse=True)
    new_verts = cyl.vertices[unique_vid]
    new_faces = inverse.reshape((-1,3))
    new_mesh = trimesh.Trimesh(vertices=new_verts, faces=new_faces, process=False)
    return new_mesh

# ---- popravljena funkcija ----
def clip_segment_to_plane(p0, p1, z_clip=0.0):
    """Ispravno odseca segment na ravni z=z_clip"""
    z0 = p0[2]; z1 = p1[2]

    if z0 >= z_clip - Z_TOL and z1 >= z_clip - Z_TOL:
        return (p0, p1)        # oba iznad
    if z0 < z_clip - Z_TOL and z1 < z_clip - Z_TOL:
        return None            # oba ispod

    # presecanje
    t = (z_clip - z0) / (z1 - z0 + 1e-16)
    t = np.clip(t, 0.0, 1.0)
    p_int = p0 + t * (p1 - p0)

    if z0 >= z_clip - Z_TOL:
        return (p0, p_int)     # gornja tačka do preseka
    else:
        return (p_int, p1)     # presečna tačka do gornje

# ---------------------------------

def main():
    print("Generišem icosahedron i subdividujem (freq = {})...".format(freq))
    V, F = make_icosahedron_unit()
    pts, faces = subdivide_and_project(V, F, freq, R)

    Ry = rotation_matrix_y(ROT_Y_DEG)
    pts_rot = (pts @ Ry.T)

    z_clip = compute_z_clip(R, CLIP_MODE, Z_CLIP_HEIGHT_MANUAL)
    print(f"Rotacija oko Y za {ROT_Y_DEG:.6f}°, CLIP_MODE={CLIP_MODE}, z_clip={z_clip:.2f} cm")

    edges = build_edges(faces)

    edge_labels = {}
    for e in edges:
        L = edge_length(pts_rot, e)
        edge_labels[e] = nearest_label(L, GIVEN)

    cnt = Counter(edge_labels[e] for e in edge_labels)
    print("Broj po tipovima (pre clip-a):", cnt)

    colors_rgb = {'A':[200,50,50], 'B':[50,200,50], 'C':[50,50,200]}
    dxf_colors = {'A':1, 'B':3, 'C':5}  # 1=red,3=green,5=blue

    cylinders = []
    cyl_colors = []
    centerlines = []
    for i, e in enumerate(edge_labels):
        u,v = e
        p0 = pts_rot[u]
        p1 = pts_rot[v]
        lab = edge_labels[e]

        cyl_piece = cylinder_between_clipped(p0, p1, strut_radius, cylinder_sections, z_clip)
        if cyl_piece is not None:
            cylinders.append(cyl_piece)
            cyl_colors.append(colors_rgb.get(lab, [150,150,150]))
        seg = clip_segment_to_plane(p0, p1, z_clip=z_clip)
        if seg is not None:
            centerlines.append((seg[0], seg[1], lab))

    if not cylinders:
        print("Nema cilindara iznad z_clip.")
        return

    combined = trimesh.util.concatenate(cylinders)

    vert_offset = 0
    verts_list = []
    faces_list = []
    colors_list = []
    for idx, cyl in enumerate(cylinders):
        nv = len(cyl.vertices)
        faces_list.append(cyl.faces + vert_offset)
        verts_list.append(cyl.vertices)
        color4 = np.tile(np.array(cyl_colors[idx] + [255], dtype=np.uint8), (nv,1))
        colors_list.append(color4)
        vert_offset += nv

    all_verts = np.vstack(verts_list)
    all_faces = np.vstack(faces_list)
    all_colors = np.vstack(colors_list)

    mesh_with_color = trimesh.Trimesh(vertices=all_verts, faces=all_faces, process=False)
    mesh_with_color.visual.vertex_colors = all_colors

    mesh_with_color.export(OUT_PLY)
    combined.export(OUT_STL)

    # --- DXF export sa bojama po tipu ---
    doc = ezdxf.new(dxfversion='R2010')
    msp = doc.modelspace()
    for (p0,p1,lab) in centerlines:
        if lab not in doc.layers:
            doc.layers.new(name=lab, dxfattribs={'color': dxf_colors.get(lab, 7)})
        msp.add_line(tuple(p0.tolist()), tuple(p1.tolist()), dxfattribs={'layer': lab})
    doc.saveas(OUT_DXF)

    print("Gotovo. Sačuvano:", OUT_PLY, OUT_STL, OUT_DXF)

if __name__ == "__main__":
    main()