Мобильная версияЖелаемая форма поверхности:
Реконструированная поверхность с помощью Делоне:
#!/usr/bin/env python3
"""
RTSTRUCT ROI triangulation (MATLAB delaunayTriangulation + freeBoundary equivalent)
Requirements:
pip install pydicom numpy scipy
Optional (for saving):
pip install trimesh pyglet
Usage:
python triangulate.py --rtstruct path/to/file.dcm --roi TECH_PTV_6600 --out ptv_mesh.ply
"""
from __future__ import annotations
import argparse
from dataclasses import dataclass
from typing import Tuple, List, Optional, Dict
import numpy as np
import pydicom
from scipy.spatial import Delaunay
class SurfaceMesh:
def __init__(self, vertices, faces):
self.vertices = vertices # (V, 3) float64
self.faces = faces # (F, 3) int64 (triangle vertex indices)
def load_dicom(path: str):
return pydicom.dcmread(path)
def build_roi_name_to_number(rs) -> Dict[str, int]:
"""
StructureSetROISequence entries have ROIName and ROINumber. This function saves the ROIName with corresponding ROINumber
"""
roi_map = {}
if not hasattr(rs, "StructureSetROISequence"):
raise ValueError("Not an RTSTRUCT: missing StructureSetROISequence.")
for roi in rs.StructureSetROISequence:
roi_map[str(roi.ROIName)] = int(roi.ROINumber)
return roi_map
def find_roi_contour_item(rs, roi_number: int):
"""
Find the ROIContourSequence item that references roi_number.
"""
if not hasattr(rs, "ROIContourSequence"):
raise ValueError("RTSTRUCT missing ROIContourSequence.")
for item in rs.ROIContourSequence:
if int(item.ReferencedROINumber) == int(roi_number):
return item
raise ValueError(f"ROIContourSequence item not found for ROINumber={roi_number}.")
def extract_roi_point_cloud_mm(rs, roi_name: str) -> np.ndarray:
"""
Collect all contour points (x,y,z) for the ROI into a single (N,3) array in mm.
This is the point cloud used by the MATLAB code as input to delaunayTriangulation.
"""
roi_map = build_roi_name_to_number(rs)
if roi_name not in roi_map:
available = "\n".join(sorted(roi_map.keys()))
raise ValueError(f'ROI "{roi_name}" not found. Available ROIs:\n{available}')
roi_number = roi_map[roi_name]
roi_contour = find_roi_contour_item(rs, roi_number)
if not hasattr(roi_contour, "ContourSequence"):
raise ValueError(f'ROI "{roi_name}" has no ContourSequence.')
pts_all = []
for c in roi_contour.ContourSequence:
data = np.asarray(c.ContourData, dtype=np.float64).reshape(-1, 3) # extract all the points of the current slice into data in format (N, 3)
pts_all.append(data) # append array into pts_all
pts = np.vstack(pts_all) # stacking all the points into a single (N, 3) array
# Remove exact duplicates (helps Delaunay stability)
pts = np.unique(np.round(pts, decimals=6), axis=0)
if pts.shape[0] < 4:
raise ValueError("Not enough unique points to build a 3D triangulation.")
return pts
def free_boundary_from_tetrahedra(delaunay: Delaunay) -> np.ndarray:
tets = delaunay.simplices # (T,4)
f0 = tets[:, [0, 1, 2]]
f1 = tets[:, [0, 1, 3]]
f2 = tets[:, [0, 2, 3]]
f3 = tets[:, [1, 2, 3]]
faces = np.vstack([f0, f1, f2, f3]) # (4T,3)
faces_sorted = np.sort(faces, axis=1)
faces_sorted = np.ascontiguousarray(faces_sorted)
unique_faces, counts = np.unique(faces_sorted, axis=0, return_counts=True)
boundary_faces = unique_faces[counts == 1]
return boundary_faces.astype(np.int64)
def triangulate_point_cloud_surface(points_mm: np.ndarray) -> SurfaceMesh:
"""
Equivalent of:
DT = delaunayTriangulation(x,y,z)
[Faces,Vertices] = freeBoundary(DT)
"""
# Delaunay tetrahedralization
dt = Delaunay(points_mm)
# Boundary triangles
faces = free_boundary_from_tetrahedra(dt)
# Vertices are just the Delaunay input points
vertices = np.asarray(dt.points, dtype=np.float64)
return SurfaceMesh(vertices=vertices, faces=faces)
def save_mesh(mesh: SurfaceMesh, out_path: str) -> None:
"""
Optional: save as PLY/OBJ/STL using trimesh if available.
"""
try:
import trimesh
except ImportError as e:
raise RuntimeError(
"Saving requires trimesh. Install with: pip install trimesh"
) from e
tm = trimesh.Trimesh(vertices=mesh.vertices, faces=mesh.faces, process=False)
#tm.show()
tm.export(out_path)
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--rtstruct", default="U:\Data\RS1.2.752.243.1.1.20260211130851705.2580.11322.dcm", help="Path to RTSTRUCT DICOM (.dcm) file")
ap.add_argument("--roi", default="TECH_PTV_6600", help='ROI name, e.g. "TECH_PTV_6600"')
ap.add_argument("--out", default="ptv_mesh.ply", help="Optional: output mesh file (ply/obj/stl)")
args = ap.parse_args()
rs = load_dicom(args.rtstruct)
pts = extract_roi_point_cloud_mm(rs, args.roi)
mesh = triangulate_point_cloud_surface(pts)
print(f"ROI: {args.roi}")
print(f"Point cloud: {pts.shape[0]} points")
print(f"Surface mesh: {mesh.vertices.shape[0]} vertices, {mesh.faces.shape[0]} triangles")
print(f"Surface mesh: {len(mesh.vertices)} vertices, {len(mesh.faces)} triangles")
if args.out:
save_mesh(mesh, args.out)
print(f"Saved mesh to: {args.out}")
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(mesh.vertices[:,0],
mesh.vertices[:,1],
mesh.vertices[:,2],
triangles=mesh.faces,
linewidth=0.2,
alpha=0.7)
plt.show()
if __name__ == "__main__":
main()
Подробнее здесь: https://stackoverflow.com/questions/798 ... ted-suface
