diff --git a/.gitignore b/.gitignore index 112ee266..b3be9d63 100644 --- a/.gitignore +++ b/.gitignore @@ -36,4 +36,53 @@ Zc*.mat *.stl ChatGPT*.html -ChatGPT*.txt \ No newline at end of file +ChatGPT*.txt + +# OS-generated files +.DS_Store + +# Compiled binaries and object files +*.o +*.a +*.so +*.dll +*.exe +*.out + +# Build output directories (example, adjust as needed) +/bin/ +/build/ +/dist/ + +# Other common temporary files +*.tmp +*.npy + +# Generated geometry +*.stl +*.step +*.stp +*.obj + +# PCB / fabrication outputs +*.gbr +*.gerber +*.drl + +# Simulation / intermediate outputs +*.vtk +*.vtu +*.npy +*.npz + +# Plot outputs +*.png +*.jpg +*.svg + +# Temporary / logs +*.log +tmp/ +output/ +build/ +dist/ \ No newline at end of file diff --git a/docs/source/api_new_modules.rst b/docs/source/api_new_modules.rst new file mode 100644 index 00000000..bb9f43be --- /dev/null +++ b/docs/source/api_new_modules.rst @@ -0,0 +1,35 @@ +New Modules API +================= + +Metrics +------- + +.. automodule:: pyCoilGen.helpers.metrics + :members: + +Plotting +-------- + +.. automodule:: pyCoilGen.plotting.plot_field + :members: + +.. automodule:: pyCoilGen.plotting.plot_wire_loops + :members: + +Wire Extraction +--------------- + +.. automodule:: pyCoilGen.sub_functions.extract_wire_paths + :members: + +Gradient Former +---------------- + +.. automodule:: pyCoilGen.sub_functions.gradient_former + :members: + +Simulation +----------- + +.. automodule:: pyCoilGen.sub_functions.simulate_gradient_coil + :members: diff --git a/docs/source/index.rst b/docs/source/index.rst index 1ae79beb..0b1dbeb0 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -34,4 +34,15 @@ pyCoilGen User Guide quick_start.md configuration.md results.md - glossary.md \ No newline at end of file + glossary.md + +What's New +========== + +This update introduces a significantly more Pythonic and modular workflow for gradient coil design and analysis: + +- Native Python-based STL generation +- Automated wire path extraction +- Integrated gradient coil simulation utilities +- New plotting utilities +- Modular helper functions for metrics \ No newline at end of file diff --git a/docs/source/migration.rst b/docs/source/migration.rst new file mode 100644 index 00000000..b1ae5a1d --- /dev/null +++ b/docs/source/migration.rst @@ -0,0 +1,15 @@ +Migration Guide +=============== + +Key Changes +----------- + +- Python-based STL generation +- Explicit wire extraction +- Integrated simulation + +Updates Required +---------------- + +- Use extract_wire_paths +- Replace geometry loading with generate_gradient_former diff --git a/docs/source/workflow.rst b/docs/source/workflow.rst new file mode 100644 index 00000000..19fcabf2 --- /dev/null +++ b/docs/source/workflow.rst @@ -0,0 +1,36 @@ +Updated Workflow for Planar Coil Design +====================================== + +The coil design pipeline has been refactored into a more modular and Python-native workflow. + +Pipeline Overview +----------------- + +1. Define geometry and parameters +2. Solve stream function optimization +3. Extract wire paths +4. Generate gradient former geometry +5. Simulate magnetic fields +6. Visualize results + +Example Workflow +---------------- + +.. code-block:: python + + from pyCoilGen.sub_functions.stream_function_optimization import optimize_stream_function + from pyCoilGen.sub_functions.extract_wire_paths import extract_wire_paths + from pyCoilGen.sub_functions.gradient_former import generate_gradient_former + from pyCoilGen.sub_functions.simulate_gradient_coil import simulate_gradient_coil + + # Step 1: Optimize stream function + sf = optimize_stream_function(...) + + # Step 2: Extract wire paths + loops = extract_wire_paths(sf, ...) + + # Step 3: Generate printable geometry + generate_gradient_former(loops, ...) + + # Step 4: Simulate field + field = simulate_gradient_coil(loops, ...) \ No newline at end of file diff --git a/examples/biplanar_gradient.py b/examples/biplanar_gradient.py new file mode 100644 index 00000000..257aeb75 --- /dev/null +++ b/examples/biplanar_gradient.py @@ -0,0 +1,181 @@ +import sys +sys.path.append('.') +from pyCoilGen.sub_functions.constants import DEBUG_BASIC, DEBUG_VERBOSE +import numpy as np +import trimesh +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +import trimesh +import logging +from pyCoilGen.pyCoilGen_release import pyCoilGen +from pyCoilGen.pyCoilGen_planar import pyCoilGen_planar +from scipy.special import sph_harm_y +# from pyCoilGen.sub_functions.utils import * +from pyCoilGen.sub_functions.stl_mesh_generation import * +from pyCoilGen.sub_functions.simulate_gradient_coil import simulate_gradient_coil +from pyCoilGen.sub_functions.gradient_former import generate_gradient_former +from pyCoilGen.helpers.metrics import print_metrics + +logging.basicConfig(level=logging.INFO) +log = logging.getLogger(__name__) + +# ========================================================== +# MAIN PROGRAM +# ========================================================== + +if __name__ == "__main__": + + """Planar MRI gradient coil design pipeline using pyCoilGen. + This script demonstrates a complete end-to-end workflow for designing, + simulating, evaluating, and exporting a biplanar MRI gradient coil. + Workflow: + 1. Create and save an STL mesh for the current-carrying planar surfaces. + 2. Validate mesh quality and visualize geometry. + 3. Configure pyCoilGen optimization and discretization parameters. + 4. Run planar coil synthesis with ``pyCoilGen_planar``. + 5. Simulate magnetic field behavior over the DSV with Magpylib helpers. + 6. Compute and print gradient performance metrics. + 7. Export coil former and wire-channel STL files for fabrication. + Attributes: + PLATE_SIZE (float): Plate width/height in meters (6 in). + PLATE_TOP (float): Z position of the top plate in meters. + PLATE_BOTTOM (float): Z position of the bottom plate in meters. + DSV_IMAGING (float): Diameter of the spherical imaging region in meters. + CNC (dict[str, float]): Fabrication and conductor parameters: + - ``diameter``: Conductor diameter in meters. + - ``current``: Current per contour in amperes. + - ``thickness``: Wire thickness / cut width in meters. + - ``width``: Wire width / normal shift length in meters. + MESH_RESOLUTION (int): General mesh-resolution control value. + Notes: + - The target field is configured through ``field_shape_function`` (e.g., "z"). + - ``target_gradient_strength`` is interpreted in T/m at 1 A drive current. + - Output artifacts (plots, debug data, exported STL files) are written to + the configured project/output directories. + ============================================================== + PLANAR MRI GRADIENT COIL DESIGN PIPELINE + ============================================================== + This example demonstrates the design of a biplanar gradient coil using pyCoilGen. The pipeline includes the following steps: + 1. Design the STL file for the current-carrying surface and save it. + 2. Check the STL file quality and visualize the geometry. + 3. Configure pyCoilGen inputs for the coil design. + 4. Run pyCoilGen to obtain the coil design and wire patterns. + 5. Simulate the coil design using Magpylib to obtain the magnetic field distribution. + 6. Evaluate the performance metrics of the gradient coil design. + 7. Write the STL files for the coil parts (plates and wire channels) for CNC machining or 3D printing. + """ + + log.info("\n========================================") + log.info("PLANAR GRADIENT COIL DESIGN PIPELINE") + log.info("========================================\n") + + # ========================================================== + # GLOBAL PARAMETERS + # ========================================================== + + PLATE_SIZE = 152.4e-3 # in meters, 6 inches + PLATE_TOP = 0.03675 + PLATE_BOTTOM = -0.03675 + DSV_IMAGING = 0.032 # in meters, diameter of the spherical imaging region + CNC = { + "diameter": 2e-3, # meters + "current": 1.0, # current per contour (amps) + "thickness": 1.6e-3, # meters, thickness of the wire (overloaded as interconnection cut width) + "width": 2e-3, # meters, width of the wire (overloaded as normal shift length) + } + MESH_RESOLUTION = 40 + + # ========================================================== + # STEP 1 - DESIGN THE STL FILE FOR CURRENT CARRYING SURFACE - SAVE STL + # ========================================================== + stl_path = "data/pyCoilGenData/Geometry_Data/Tenacity_circular_1.stl" + radius = 0.07 + nr = 30 # number of radial divisions + nt = 90 # 180 for 2 degree resolution, 90 for 4 degree resolution + inner_radius = radius / nr + z_positions = [-0.03675, 0.03675] + create_stl_mesh(radius,nr,nt,z_positions,stl_path) + + #========================================================= + # STEP 2 - CHECK THE STL FILE QUALITY AND VISUALIZE THE GEOMETRY + #======================================================== + mesh = trimesh.load(stl_path) + check_mesh_quality(mesh, dsv_radius=DSV_IMAGING * 0.5) + + #========================================================= + # STEP 3 - CONFIGURE PYCOILGEN INPUTS + #========================================================= + arg_dict = { + 'field_shape_function': 'z', # definition of the target field + 'coil_mesh_file': stl_path.split('/')[-1],# assumes the STL file is in the Geometry_Data folder + # 'target_mesh_file': None, + 'secondary_target_mesh_file': 'none', + # 'secondary_target_weight': 0.5, + 'target_region_radius': DSV_IMAGING * 0.5, # in meter + # 'target_region_resolution': 10, # MATLAB 10 is the default + 'target_region_resolution': 10, # number of points along each axis in the target region; higher values lead to better optimization results but longer runtime + 'use_only_target_mesh_verts': False, + 'sf_source_file': 'none', + # the number of potential steps that determines the later number of windings (Stream function discretization) + 'levels': 20, # was 12 + # a potential offset value for the minimal and maximal contour potential ; must be between 0 and 1 + 'pot_offset_factor': 0.35, + 'surface_is_cylinder_flag': False, + # the width for the interconnections are interconnected; in meter + 'interconnection_cut_width': CNC['thickness'], # overloaded as wire_depth + # the length for which overlapping return paths will be shifted along the surface normals; in meter + 'normal_shift_length': CNC['width'], # overloaded as wire_spacing = 2mm for AWG14 and 0.3mm for CNC + 'iteration_num_mesh_refinement': 1, # the number of refinements for the mesh; + 'set_roi_into_mesh_center': True, + 'force_cut_selection': ['high'], + # Specify one of the three ways the level sets are calculated: "primary","combined", or "independent" + 'skip_postprocessing': False, + 'skip_inductance_calculation': False, + 'tikhonov_reg_factor': 200, # Tikhonov regularization factor for the SF optimization + 'output_directory': 'images', # [Current directory] + 'project_name': 'biplanar_xgradient_Tenacity', + 'persistence_dir': 'debug', + 'debug': DEBUG_VERBOSE, + 'target_field_weighting': True, + 'minimize_method_parameters': "{'tol': 1e-9}", + 'minimize_method_options': "{'disp': True, 'maxiter' : 5}", + 'target_gradient_strength': 1, # in T/m, this is the target gradient strength at the center of the DSV; the optimization will try to achieve this strength with 1 A of current + 'skip_postprocessing': True, + 'level_set_method': 'primary', + 'sf_opt_method': 'Iterative', + # 'sf_dest_file': 'test_grad_design_preopt', + # 'sf_source_file': 'test_grad_design_preopt', + # 'sf_source_file': 'none', + } + + # ========================================================== + # STEP 4 - RUN PYCOILGEN TO OBTAIN THE COIL DESIGN AND WIRE PATTERNS + # ========================================================== + log.info("Running pyCoilGen planar...") + coil_parts = pyCoilGen_planar(log, arg_dict) + log.info("pyCoilGen finished") + + #========================================================= + # STEP 5 - SIMULATE THE COIL DESIGN USING MAGPYLIB TO OBTAIN THE MAGNETIC FIELD DISTRIBUTION + #========================================================= + # simulation contains - B, Bx, By, Bz, points (where the field is evaluated), and the coil geometry + simulation = simulate_gradient_coil(coil_parts, DSV_IMAGING, wire=CNC, display_field=True) + + # ======================================================== + # STEP 6 - EVALUATE THE PERFORMANCE METRICS OF THE GRADIENT COIL DESIGN + # ======================================================== + log.info("Evaluating performance metrics...") + print_metrics(simulation['Bz'], coords = simulation['points'], axis = arg_dict['field_shape_function']) + + #========================================================== + # STEP 7 - WRITE THE STL FILES FOR THE COIL PARTS (PLATES AND WIRE CHANNELS) FOR CNC MACHINING + #========================================================== + output_prefix = f"./images/{arg_dict['project_name']}_{arg_dict['field_shape_function']}" + generate_gradient_former( + coil_parts, + output_prefix=output_prefix) + + + log.info("\n========================================") + log.info("PLANAR GRADIENT COIL DESIGN PIPELINE COMPLETED") + log.info("========================================\n") \ No newline at end of file diff --git a/pyCoilGen/helpers/metrics.py b/pyCoilGen/helpers/metrics.py new file mode 100644 index 00000000..641d4a8f --- /dev/null +++ b/pyCoilGen/helpers/metrics.py @@ -0,0 +1,93 @@ +import numpy as np +import matplotlib.pyplot as plt + + +# ========================================================== +# PRINT PERFORMANCE METRICS +# ========================================================== + +def print_metrics(Bz, coords, axis = 'x'): + """ + Print and visualize magnetic gradient field metrics. + This function creates a 3D scatter plot of the magnetic field and computes + gradient efficiency and linearity error metrics along the specified axis. + Args: + Bz (numpy.ndarray): Magnetic field values in Tesla. 1D array of field + measurements at each coordinate point. + coords (numpy.ndarray): Coordinate positions. 2D array of shape (N, 3) + containing x, y, z coordinates in meters for each field measurement. + axis (str, optional): Axis along which to evaluate gradient metrics. + Must be one of 'x', 'y', or 'z'. Defaults to 'x'. + Returns: + None + Raises: + ValueError: If axis is not one of 'x', 'y', or 'z'. + Prints: + - Gradient efficiency (mT/m/A): Linear fit slope along specified axis. + - Max linearity error within DSV (%): Maximum deviation from linear fit + expressed as a percentage of the maximum fitted value. + Note: + The function assumes the measurement points within the designated + sampling volume (DSV) are those within ±1 mm of the respective axis. + """ + + + x = coords[:, 0] + y = coords[:, 1] + z = coords[:, 2] + + fig = plt.figure(figsize=(8, 6)) + ax = fig.add_subplot(111, projection='3d') + sc = ax.scatter(x, y, z, c=Bz * 1e3, cmap='jet') + plt.colorbar(sc, label='Gradient Field (mT)') + ax.set_title('Optimized: ' + axis + ' Magnetic Gradient Field') + ax.set_xlabel('X (m)') + ax.set_ylabel('Y (m)') + ax.set_zlabel('Z (m)') + plt.show() + + # Evaluate along X-axis + if axis == 'x': + axis_mask = (np.abs(y) < 1e-3) & (np.abs(z) < 1e-3) + x_axis = x[axis_mask] + g_axis = Bz[axis_mask] + + coeffs = np.polyfit(x_axis, g_axis, 1) + slope = coeffs[0] + + print("\n===== X GRADIENT PERFORMANCE =====") + print(f"Gradient efficiency (mT/m/A): {slope*1e3:.3f}") + + linear_fit = np.polyval(coeffs, x_axis) + error = (g_axis - linear_fit) / np.max(np.abs(linear_fit)) * 100 + print(f"Max linearity error within DSV (%): {np.max(np.abs(error)):.2f}") + + elif axis == 'y': + axis_mask = (np.abs(x) < 1e-3) & (np.abs(z) < 1e-3) + y_axis = y[axis_mask] + g_axis = Bz[axis_mask] + + coeffs = np.polyfit(y_axis, g_axis, 1) + slope = coeffs[0] + + print("\n===== Y GRADIENT PERFORMANCE =====") + print(f"Gradient efficiency (mT/m/A): {slope*1e3:.3f}") + + linear_fit = np.polyval(coeffs, y_axis) + error = (g_axis - linear_fit) / np.max(np.abs(linear_fit)) * 100 + print(f"Max linearity error within DSV (%): {np.max(np.abs(error)):.2f}") + + elif axis == 'z': + axis_mask = (np.abs(x) < 1e-3) & (np.abs(y) < 1e-3) + z_axis = z[axis_mask] + g_axis = Bz[axis_mask] + + coeffs = np.polyfit(z_axis, g_axis, 1) + slope = coeffs[0] + + print("\n===== Z GRADIENT PERFORMANCE =====") + print(f"Gradient efficiency (mT/m/A): {slope*1e3:.3f}") + + linear_fit = np.polyval(coeffs, z_axis) + error = (g_axis - linear_fit) / np.max(np.abs(linear_fit)) * 100 + print(f"Max linearity error within DSV (%): {np.max(np.abs(error)):.2f}") \ No newline at end of file diff --git a/pyCoilGen/plotting/plot_field.py b/pyCoilGen/plotting/plot_field.py new file mode 100644 index 00000000..08c08640 --- /dev/null +++ b/pyCoilGen/plotting/plot_field.py @@ -0,0 +1,36 @@ +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + + +def plot_field(Bx, By, Bz, points): + """Visualize magnetic field components in 3D space. + Creates a figure with three subplots displaying the x, y, and z components + of a magnetic field. Each component is visualized as a 3D scatter plot with + color mapping representing the field magnitude at each point. + Args: + Bx (array-like): Magnetic field component in the x-direction. + By (array-like): Magnetic field component in the y-direction. + Bz (array-like): Magnetic field component in the z-direction. + points (ndarray): Array of shape (N, 3) containing the 3D coordinates + of points where the field is sampled, with columns representing + x, y, and z coordinates respectively. + Returns: + None + Raises: + None + Note: + The function displays the plot using matplotlib's interactive backend. + All three field components must have the same length as the number of rows + in the points array. + """ + + fig = plt.figure(figsize=(15,5)) + for i, (comp, title) in enumerate(zip([Bx, By, Bz], ["Bx","By","Bz"])): + ax = fig.add_subplot(1, 3, i+1, projection='3d') + sc = ax.scatter(points[:,0], points[:,1], points[:,2], c=comp, cmap='jet') + ax.set_title(title) + fig.colorbar(sc, ax=ax, shrink=0.6) + fig.suptitle("Simulated Gradient Field") + plt.tight_layout() + plt.show() \ No newline at end of file diff --git a/pyCoilGen/plotting/plot_wire_loops.py b/pyCoilGen/plotting/plot_wire_loops.py new file mode 100644 index 00000000..af3cca30 --- /dev/null +++ b/pyCoilGen/plotting/plot_wire_loops.py @@ -0,0 +1,268 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle +import trimesh +from scipy.interpolate import splprep, splev +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle, Circle, FancyArrow +import pyvista as pv + +# ---------------------------- +# debug plot +# ---------------------------- +def plot_wire_loops(wire_loops): + """ + Plot wire loops for gradient coils. + Visualizes the paths of wire loops used in gradient coil design, + with different colors representing positive and negative loops. + Args: + wire_loops (list): A list of dictionaries containing wire loop data. + Each dictionary should have: + - "path" (numpy.ndarray): An (n, 2) array of (x, y) coordinates + representing the wire loop path. + - "sign" (int or float): The polarity of the loop. Positive values + are plotted in red, non-positive values in blue. + Returns: + None + Raises: + None + """ + + fig, ax = plt.subplots(figsize=(6,6)) + for loop in wire_loops: + p = loop["path"] + ax.plot(p[:,0], p[:,1], 'r' if loop["sign"]>0 else 'b') + ax.set_aspect("equal") + ax.set_title("Gradient Coil Wire Paths") + plt.show() + +def plot_wire_loops_tube(wire_loops, wire_width): + """ + Visualize wire loops as 3D tubes in a plotter. + Renders wire loops as tubular meshes with color coding based on their sign. + Positive sign loops are displayed in red, while negative sign loops are + displayed in blue. The visualization includes axes and grid for reference. + Args: + wire_loops (list): A list of dictionaries containing wire loop data. + Each dictionary must have: + - "path" (array-like): The 3D coordinates of points defining the loop path. + - "sign" (float): The sign value determining the tube color (positive + for red, non-positive for blue). + wire_width (float): The diameter of the wire tubes in plotting units. + The tube radius is calculated as wire_width/2. + Returns: + None + Raises: + None + Example: + >>> wire_loops = [ + ... {"path": [[0, 0, 0], [1, 0, 0], [1, 1, 0]], "sign": 1}, + ... {"path": [[0, 0, 1], [1, 0, 1], [1, 1, 1]], "sign": -1} + ... ] + >>> plot_wire_loops_tube(wire_loops, wire_width=0.1) + """ + + + plotter = pv.Plotter() + + for loop in wire_loops: + + poly = pv.lines_from_points(loop["path"]) + + tube = poly.tube(radius=wire_width/2) + + plotter.add_mesh( + tube, + color="red" if loop["sign"] > 0 else "blue" + ) + + plotter.add_axes() + plotter.show_grid() + plotter.show() + + +def plot_gerber_paths(gerber_paths, plate_radius_mm, wire_width_mm, part_index): + """Visualize Gerber paths with wire geometry on a circular plate. + Plots wire loops extracted from Gerber files, displaying both the path + centerlines and the actual wire cross-sections with specified width. + Positive and negative paths are distinguished by color (red and blue). + Args: + gerber_paths (list): List of tuples containing (path, sign) pairs where + path is an Nx2 numpy array of (x, y) coordinates and sign indicates + the polarity (+1 or -1) of the path. + plate_radius_mm (float): Radius of the circular plate boundary in + millimeters. Displayed as a dashed gray circle. + wire_width_mm (float): Width of the wire in millimeters. Used to + render rectangular cross-sections along each path segment. + part_index (int): Index or identifier of the plate being plotted. + Displayed in the figure title. + Returns: + None + Note: + - Red paths represent positive polarity (sign > 0) + - Blue paths represent negative polarity (sign <= 0) + - Wire cross-sections are rendered as rotated rectangles aligned + with path segments + - The plot uses equal aspect ratio for accurate geometric representation + """ + + plt.figure(figsize=(6,6)) + ax = plt.gca() + + for path, sign in gerber_paths: + + ax.plot( + path[:,0], + path[:,1], + color="red" if sign>0 else "blue", + linewidth=1.5 + ) + + for i in range(len(path)-1): + + p0, p1 = path[i], path[i+1] + + dx, dy = p1[0]-p0[0], p1[1]-p0[1] + + length = np.sqrt(dx**2 + dy**2) + + angle = np.degrees(np.arctan2(dy, dx)) + + + + rect = Rectangle( + (p0[0]-wire_width_mm/2, p0[1]-wire_width_mm/2), + length, + wire_width_mm, + angle=angle, + color="red" if sign>0 else "blue", + alpha=0.2 + ) + + ax.add_patch(rect) + + circle = plt.Circle( + (0,0), + plate_radius_mm, + fill=False, + linestyle='--', + color='gray' + ) + + ax.add_patch(circle) + + ax.set_aspect('equal') + ax.set_xlabel("X [mm]") + ax.set_ylabel("Y [mm]") + + plt.title( + f"Plate {part_index} Wire Loops with {wire_width_mm:.1f} mm Width" + ) + + plt.show() + +def plot_stl_patch(gerber_paths, plate_radius_mm, wire_width_mm, part_index): + """ + Plot wire loops with terminals and strain relief on a circular plate. + Visualizes the layout of wire loops on a circular plate, including the plate + boundary, wire paths, grooves, terminal holes, and strain relief indicators. + Wire paths are color-coded by sign (red for positive, blue for negative). + Args: + gerber_paths (list of tuples): List of (path, sign) tuples where path is + a numpy array of shape (n, 2) containing 2D coordinates and sign is + a numeric value indicating the loop polarity (positive or negative). + plate_radius_mm (float): Radius of the circular plate in millimeters. + wire_width_mm (float): Width of the wire in millimeters, used for + calculating groove dimensions and terminal hole sizes. + part_index (int): Index or identifier of the plate for display in the + plot title. + Returns: + None: Displays the plot using matplotlib. + Raises: + None + Note: + - Red wire paths indicate positive polarity, blue indicate negative. + - Green semi-transparent circles mark the start and end terminal holes. + - Orange arrows indicate strain relief directions. + - Groove widths are visualized as semi-transparent rectangles along + the wire paths. + """ + + plt.figure(figsize=(8,8)) + ax = plt.gca() + # Draw plate boundary + circle = Circle( + (0,0), + plate_radius_mm, + fill=False, + linestyle='--', + color='gray', + linewidth=2 + ) + ax.add_patch(circle) + + # Draw wire loops + for loop_index, (path, sign) in enumerate(gerber_paths): + ax.plot( + path[:,0], + path[:,1], + color="red" if sign>0 else "blue", + linewidth=1.5, + label=f"Loop {loop_index+1}" + ) + + # Draw grooves as semi-transparent rectangles + for i in range(len(path)-1): + p0, p1 = path[i], path[i+1] + dx, dy = p1-p0 + length = np.sqrt(dx**2 + dy**2) + angle = np.degrees(np.arctan2(dy, dx)) + rect = Rectangle( + (p0[0]-wire_width_mm/2, p0[1]-wire_width_mm/2), + length, + wire_width_mm, + angle=angle, + color="red" if sign>0 else "blue", + alpha=0.2 + ) + ax.add_patch(rect) + + # Draw start and end terminal holes + def tangent(p0, p1): + v = p1 - p0 + n = np.linalg.norm(v) + return v/n if n>0 else np.array([1.0,0.0]) + def perp(v): return np.array([-v[1], v[0]]) + + t_start = tangent(path[0], path[1]) + t_end = tangent(path[-2], path[-1]) + start_pt = path[0] + perp(t_start)*wire_width_mm*1.5 + end_pt = path[-1] + perp(t_end)*wire_width_mm*1.5 + + ax.add_patch(Circle(start_pt, radius=wire_width_mm*0.8, color='green', alpha=0.5)) + ax.add_patch(Circle(end_pt, radius=wire_width_mm*0.8, color='green', alpha=0.5)) + + # Draw strain relief arrows + strain_len = wire_width_mm*4 + ax.add_patch(FancyArrow( + start_pt[0], start_pt[1], + t_start[0]*strain_len, t_start[1]*strain_len, + width=0.3, color='orange', alpha=0.6 + )) + ax.add_patch(FancyArrow( + end_pt[0], end_pt[1], + t_end[0]*strain_len, t_end[1]*strain_len, + width=0.3, color='orange', alpha=0.6 + )) + + # Draw loop numbers + ax.text(path[0][0], path[0][1], str(loop_index+1), + fontsize=10, color='purple', weight='bold') + + ax.set_aspect('equal') + ax.set_xlabel("X [mm]") + ax.set_ylabel("Y [mm]") + plt.title(f"Plate {part_index} Wire Loops with Terminals and Strain Relief") + plt.grid(True) + plt.show() \ No newline at end of file diff --git a/pyCoilGen/pyCoilGen_planar.py b/pyCoilGen/pyCoilGen_planar.py new file mode 100644 index 00000000..b9242ee2 --- /dev/null +++ b/pyCoilGen/pyCoilGen_planar.py @@ -0,0 +1,367 @@ +# Logging +import logging + +# System imports +import numpy as np +from os import makedirs + +# Logging +import logging + + +# Local imports +from .sub_functions.constants import * +from .sub_functions.data_structures import CoilSolution +import matplotlib.pyplot as plt +import trimesh +import matplotlib.tri as mtri + +# For visualisation +from .helpers.visualisation import visualize_vertex_connections, visualize_compare_contours + +# For timing +from .helpers.timing import Timing + +# For saving Pickle files +from .helpers.persistence import save, save_preoptimised_data + +# From original project +from .sub_functions.read_mesh import read_mesh +from .sub_functions.parse_input import parse_input, create_input +from .sub_functions.split_disconnected_mesh import split_disconnected_mesh +from .sub_functions.refine_mesh import refine_mesh_delegated as refine_mesh +# from .sub_functions.refine_mesh import refine_mesh # Broken +from .sub_functions.parameterize_mesh import parameterize_mesh +from .sub_functions.define_target_field import define_target_field +# from .sub_functions.temp_evaluation import temp_evaluation +from .sub_functions.calculate_one_ring_by_mesh import calculate_one_ring_by_mesh +from .sub_functions.calculate_basis_functions import calculate_basis_functions +from .sub_functions.calculate_sensitivity_matrix import calculate_sensitivity_matrix +from .sub_functions.calculate_gradient_sensitivity_matrix import calculate_gradient_sensitivity_matrix +from .sub_functions.calculate_resistance_matrix import calculate_resistance_matrix +from .sub_functions.stream_function_optimization import stream_function_optimization +from .sub_functions.calc_potential_levels import calc_potential_levels +from .sub_functions.calc_contours_by_triangular_potential_cuts import calc_contours_by_triangular_potential_cuts +from .sub_functions.process_raw_loops import process_raw_loops +from .sub_functions.find_minimal_contour_distance import find_minimal_contour_distance +from .sub_functions.topological_loop_grouping import topological_loop_grouping +from .sub_functions.calculate_group_centers import calculate_group_centers +from .sub_functions.interconnect_within_groups import interconnect_within_groups +from .sub_functions.interconnect_among_groups import interconnect_among_groups +from .sub_functions.shift_return_paths import shift_return_paths +from .sub_functions.generate_cylindrical_pcb_print import generate_cylindrical_pcb_print +from .sub_functions.create_sweep_along_surface import create_sweep_along_surface +from .sub_functions.calculate_inductance_by_coil_layout import calculate_inductance_by_coil_layout +from .sub_functions.load_preoptimized_data import load_preoptimized_data +from .sub_functions.evaluate_field_errors import evaluate_field_errors +from .sub_functions.calculate_gradient import calculate_gradient +from .sub_functions.export_data import export_data, check_exporter_help +from .sub_functions.extract_wire_paths import extract_wire_paths + + +# Set up logging +log = logging.getLogger(__name__) + + +def pyCoilGen_planar(log, input_args=None): + """ + Generate optimized planar coil layouts using stream function optimization. + This module implements the pyCoilGen algorithm for designing planar gradient coils. + It reads a mesh geometry, optimizes a stream function to match a target magnetic field, + and extracts wire patterns for coil fabrication. + The workflow consists of the following main steps: + 1. Read and preprocess the input mesh (refinement, parameterization) + 2. Define target magnetic field specifications + 3. Calculate basis functions and sensitivity matrices + 4. Optimize stream function using constrained optimization + 5. Discretize the stream function into wire contours + 6. Extract and interconnect wire paths + Args: + log (logging.Logger): Logger instance for recording execution progress and errors. + input_args (dict, optional): Dictionary or parsed input arguments containing: + - project_name (str): Name of the coil project + - debug (int): Debug verbosity level + - persistence_dir (str): Directory for saving intermediate results + - output_directory (str): Directory for saving output files + - sf_source_file (str): Path to preoptimized stream function file ('none' to skip) + - sf_dest_file (str): Path to save preoptimized data ('none' to skip) + - target_field_weighting (bool): Whether to apply radial weighting to target field + - normal_shift_length (float): Wire width parameter + - interconnection_cut_width (float): Wire thickness parameter + Defaults to None, which triggers interactive input parsing. + Returns: + list or CoilSolution: List of CoilPart objects containing optimized coil geometry + and wire patterns. Returns None if mesh loading fails. If exporter help is + requested, returns empty CoilSolution object. + Raises: + Exception: Re-raises any exceptions that occur during optimization with error + context indicating the last successful runpoint. + Notes: + - Adapted from original pyCoilGen with customization for planar gradient coils + - External dependencies: NumPy, Trimesh, Matplotlib + - Timing information is logged for performance profiling + - Intermediate results can be saved/loaded for iterative design + """ + + # Create optimized coil finished coil layout + # Author: Philipp Amrein, University Freiburg, Medical Center, Radiology, Medical Physics + # 5.10.2021 + + # The following external functions were used in modified form: + # intreparc@John D'Errico (2010), @matlabcentral/fileexchange + # The non-cylindrical parameterization is taken from "matlabmesh @ Ryan Schmidt rms@dgp.toronto.edu" + # based on desbrun et al (2002), "Intrinsic Parameterizations of {Surface} Meshes", NS (2021). + # Curve intersections (https://www.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections), + # MATLAB Central File Exchange. + + + + timer = Timing() + timer.start() + + # Parse the input variables + if type(input_args) is dict: + try: + if input_args['debug'] >= DEBUG_VERBOSE: + log.debug(" - converting input dict to input type.") + except KeyError: + pass + input_parser, input_args = create_input(input_args) + elif input_args is None: + input_parser, input_args = parse_input(input_args) + else: + input_args = input_args + + set_level(input_args.debug) + + project_name = f'{input_args.project_name}' + persistence_dir = input_args.persistence_dir + image_dir = input_args.output_directory + + # Create directories if they do not exist + makedirs(persistence_dir, exist_ok=True) + makedirs(image_dir, exist_ok=True) + + # Print the input variables + # DEBUG + if get_level() >= DEBUG_VERBOSE: + log.debug('Parse inputs: %s', input_args) + + + solution = CoilSolution() + solution.input_args = input_args + + if check_exporter_help(input_args): + return solution + + try: + runpoint_tag = 'test' + + if input_args.sf_source_file == 'none': + # Read the input mesh + print('Load geometry:') + coil_mesh, target_mesh, secondary_target_mesh = read_mesh(input_args) # 01 + + if coil_mesh is None: + log.info("No coil mesh, exiting.") + timer.stop() + return None + + if get_level() >= DEBUG_VERBOSE: + log.debug(" -- vertices shape: %s", coil_mesh.get_vertices().shape) # (264,3) + log.debug(" -- faces shape: %s", coil_mesh.get_faces().shape) # (480,3) + + if get_level() >= DEBUG_VERBOSE: + log.debug(" coil_mesh.vertex_faces: %s", coil_mesh.trimesh_obj.vertex_faces[0:10]) + + if get_level() > DEBUG_VERBOSE: + coil_mesh.display() + + # Split the mesh and the stream function into disconnected pieces + print('Split the mesh and the stream function into disconnected pieces.') + timer.start() + coil_parts = split_disconnected_mesh(coil_mesh) # 00 + timer.stop() + solution.coil_parts = coil_parts + runpoint_tag = '00' + + # Upsample the mesh density by subdivision + print('Upsample the mesh by subdivision:') + timer.start() + coil_parts = refine_mesh(coil_parts, input_args) # 01 + timer.stop() + runpoint_tag = '01' + + # Parameterize the mesh + print('Parameterize the mesh:') + timer.start() + coil_parts = parameterize_mesh(coil_parts, input_args) # 02 + timer.stop() + runpoint_tag = '02' + + if get_level() >= DEBUG_VERBOSE: + # Export data + print('Exporting initial data:') + timer.start() + export_data(solution) + timer.stop() + + # Define the target field + print('Define the target field:') + timer.start() + target_field, is_suppressed_point = define_target_field( + coil_parts, target_mesh, secondary_target_mesh, input_args) + timer.stop() + # Can set target field weights here + if input_args.target_field_weighting == True: + coords = target_field.coords + x = coords[0, :] + y = coords[1, :] + z = coords[2, :] + r = np.sqrt(x**2 + y**2 + z**2) + R = np.max(r) + # weights = 1 / (1 + (r/0.01)**4) + weights = (r / R)**4 + target_field.weights = weights + + solution.target_field = target_field + solution.is_suppressed_point = is_suppressed_point + runpoint_tag = '02b' + + if get_level() >= DEBUG_VERBOSE: + log.debug(" -- target_field.b shape: %s", target_field.b.shape) # (3, 257) + log.debug(" -- target_field.coords shape: %s", target_field.coords.shape) # (3, 257) + log.debug(" -- target_field.weights shape: %s", target_field.weights.shape) # (257,) + + # Evaluate the temp data; check whether precalculated values can be used from previous iterations + # print('Evaluate the temp data:') + # input_args = temp_evaluation(solution, input_args, target_field) + + # Find indices of mesh nodes for one ring basis functions + print('Calculate mesh one ring:') + timer.start() + coil_parts = calculate_one_ring_by_mesh(coil_parts) # 03 + timer.stop() + runpoint_tag = '03' + + # Create the basis function container which represents the current density + print('Create the basis function container which represents the current density:') + timer.start() + coil_parts = calculate_basis_functions(coil_parts) # 04 + timer.stop() + runpoint_tag = '04' + + # Calculate the sensitivity matrix Cn + print('Calculate the sensitivity matrix:') + timer.start() + coil_parts = calculate_sensitivity_matrix(coil_parts, target_field, input_args) # 05 + timer.stop() + runpoint_tag = '05' + # Print the condition number of the sensitivity matrix for the first coil part + if get_level() >= DEBUG_VERBOSE: + for part_ind in range(len(coil_parts)): + coil_part = coil_parts[part_ind] + sensitivity_matrix = coil_part.sensitivity_matrix + cond_number = np.linalg.cond(sensitivity_matrix.reshape(3, -1)) + log.info("Condition number of sensitivity matrix for part %d: %e", part_ind, cond_number) + + + + + # Calculate the gradient sensitivity matrix Gn + print('Calculate the gradient sensitivity matrix:') + timer.start() + coil_parts = calculate_gradient_sensitivity_matrix(coil_parts, target_field, input_args) # 06 + timer.stop() + runpoint_tag = '06' + + # Calculate the resistance matrix Rmn + print('Calculate the resistance matrix:') + timer.start() + coil_parts = calculate_resistance_matrix(coil_parts, input_args) # 07 + timer.stop() + runpoint_tag = '07' + + # Optimize the stream function toward target field and further constraints + print('Optimize the stream function toward target field and secondary constraints:') + timer.start() + coil_parts, combined_mesh, sf_b_field = stream_function_optimization(coil_parts, target_field, input_args) + timer.stop() + solution.combined_mesh = combined_mesh + solution.sf_b_field = sf_b_field + runpoint_tag = '08' + + if input_args.sf_dest_file != 'none': + print('Persist pre-optimised data:') + save_preoptimised_data(solution) + + else: + # Load the preoptimised data + print('Load pre-optimised data:') + timer.start() + solution = load_preoptimized_data(input_args) + timer.stop() + coil_parts = solution.coil_parts + combined_mesh = solution.combined_mesh + target_field = solution.target_field + + + # Customized workflow for planar gradient design from here on, which can be adapted by the user. The following steps are based on the original workflow of pyCoilGen, but can be modified as needed. + if get_level() >= DEBUG_VERBOSE: + psis = coil_parts[0].stream_function + plt.figure() + plt.hist(psis, bins=50) + plt.title("Stream function values distribution") + plt.xlabel("Stream function value") + plt.ylabel("Frequency") + plt.show() + + + + # Calculate the potential levels for the discretization + print('Calculate the potential levels for the discretization:') + timer.start() + coil_parts, primary_surface_ind = calc_potential_levels(coil_parts, combined_mesh, input_args) # 09 + timer.stop() + solution.primary_surface_ind = primary_surface_ind + runpoint_tag = '09' + + if get_level() >= DEBUG_VERBOSE: + for part_ind in range(len(coil_parts)): + coil_part = coil_parts[part_ind] + log.info("Part %d: potential levels: %s", part_ind, coil_part.potential_level_list) + mesh = coil_parts[part_ind].coil_mesh # Assuming you have one coil part + verts = mesh.v # Nx3 + faces = mesh.f # Mx3 + psis = coil_parts[part_ind].stream_function + # Compute triangle centers + triang = mtri.Triangulation(verts[:,0], verts[:,1], faces) + + plt.figure(figsize=(6,6)) + plt.tricontourf(triang, psis, levels= coil_parts[part_ind].potential_level_list, cmap='jet') + plt.colorbar(label="Stream function ψ") + plt.axis('equal') + plt.title("Streamfunction on coil surface") + plt.show() + + # Extract wire patterns + + for part in coil_parts: + log.info("Extracting wire paths for part with %d vertices and %d faces", part.coil_mesh.get_vertices().shape[0], part.coil_mesh.get_faces().shape[0]) + part.wire_loops = extract_wire_paths(part.coil_mesh.get_vertices(), part.coil_mesh.get_faces(), part.stream_function, + part.potential_level_list, wire_width = input_args.normal_shift_length, wire_thickness=input_args.interconnection_cut_width, display_debug_plots = True) + + + + except Exception as e: + log.error("An error occurred at runpoint %s: %s", runpoint_tag, str(e)) + raise e + + + return coil_parts + + + + + + diff --git a/pyCoilGen/sub_functions/coil_verification_metrics.py b/pyCoilGen/sub_functions/coil_verification_metrics.py new file mode 100644 index 00000000..bd5e73e1 --- /dev/null +++ b/pyCoilGen/sub_functions/coil_verification_metrics.py @@ -0,0 +1,227 @@ +############################################################ +# STEP 11: VERIFY COIL FIELD AND PERFORMANCE +############################################################ + +import numpy as np +import matplotlib.pyplot as plt +import magpylib as magpy + + +############################################################ +# 11.1 Convert loops → current wires +############################################################ + +def loops_to_magpy(loops, current=1.0): + + wires = [] + + for loop in loops: + + # Ensure closed loop + if not np.allclose(loop[0], loop[-1]): + loop = np.vstack([loop, loop[0]]) + + wire = magpy.current.Line( + current=current, + vertices=loop + ) + + wires.append(wire) + + return wires + + +############################################################ +# 11.2 Create sampling points inside DSV +############################################################ + +def create_dsv_points(radius=0.016, resolution=11): + + xs = np.linspace(-radius, radius, resolution) + ys = np.linspace(-radius, radius, resolution) + zs = np.linspace(-radius, radius, resolution) + + pts = [] + + for x in xs: + for y in ys: + for z in zs: + + if x*x + y*y + z*z <= radius*radius: + pts.append([x,y,z]) + + return np.array(pts) + + +############################################################ +# 11.3 Compute magnetic field from wires +############################################################ + +def compute_field(wires, points): + + B = magpy.getB(wires, points) + + return B + + +############################################################ +# 11.4 Fit gradient tensor +############################################################ + +def fit_gradient_tensor(points, B): + + x = points[:,0] + y = points[:,1] + z = points[:,2] + + Bz = B[:,2] + + A = np.column_stack([ + np.ones(len(points)), + x, + y, + z + ]) + + coeff, *_ = np.linalg.lstsq(A, Bz, rcond=None) + + B0 = coeff[0] + Gx = coeff[1] + Gy = coeff[2] + Gz = coeff[3] + + return B0, Gx, Gy, Gz + + +############################################################ +# 11.5 Compute linearity error +############################################################ + +def compute_linearity(points, Bz, Gx): + + x = points[:,0] + + ideal = Gx * x + + error = Bz - ideal + + linearity = np.max(np.abs(error)) / np.max(np.abs(ideal)) + + return linearity * 100 + + +############################################################ +# 11.6 Visualize gradient behavior +############################################################ + +def visualize_gradient(points, B): + + x = points[:,0] + Bz = B[:,2] + + plt.figure() + + plt.scatter(x, Bz, s=10) + + plt.xlabel("x (m)") + plt.ylabel("Bz (T)") + plt.title("Generated Gradient Field") + + plt.grid() + + plt.show() + + +############################################################ +# 11.7 Visualize 2D slice of field +############################################################ + +def visualize_field_slice(points, B): + + x = points[:,0] + z = points[:,2] + Bz = B[:,2] + + plt.figure() + + plt.scatter(x*1000, z*1000, c=Bz, s=40) + + plt.xlabel("x (mm)") + plt.ylabel("z (mm)") + + plt.title("Bz field distribution") + + plt.colorbar(label="Tesla") + + plt.gca().set_aspect('equal') + + plt.show() + + +############################################################ +# 11.8 Full evaluation routine +############################################################ + +def evaluate_coil_performance(loops): + + print("\nSTEP 11: Evaluating coil field performance\n") + + print("Creating current wires...") + + wires = loops_to_magpy(loops) + + print("Generating DSV sampling points...") + + points = create_dsv_points() + + print("Total sampling points:", len(points)) + + print("Computing magnetic field...") + + B = compute_field(wires, points) + + print("Fitting gradient tensor...") + + B0, Gx, Gy, Gz = fit_gradient_tensor(points, B) + + Bz = B[:,2] + + linearity = compute_linearity(points, Bz, Gx) + + print("\n========= COIL PERFORMANCE =========") + + print("Uniform field offset B0:", B0, "T") + + print("Gradient components:") + + print(" Gx:", Gx, "T/m/A") + print(" Gy:", Gy, "T/m/A") + print(" Gz:", Gz, "T/m/A") + + print("\nGradient efficiency:") + + print(" ", Gx*1000, "mT/m/A") + + print("\nLinearity error inside DSV:") + + print(" ", linearity, "%") + + print("====================================\n") + + visualize_gradient(points, B) + + visualize_field_slice(points, B) + + performance = { + "B0": B0, + "Gx": Gx, + "Gy": Gy, + "Gz": Gz, + "linearity_error_percent": linearity + } + + return { + "performance": performance + } + + diff --git a/pyCoilGen/sub_functions/extract_wire_paths.py b/pyCoilGen/sub_functions/extract_wire_paths.py new file mode 100644 index 00000000..71211bc8 --- /dev/null +++ b/pyCoilGen/sub_functions/extract_wire_paths.py @@ -0,0 +1,268 @@ + +import numpy as np +import matplotlib.tri as mtri +import matplotlib.pyplot as plt +from scipy.interpolate import splprep, splev +from scipy.spatial import cKDTree +import pyvista as pv +from pyCoilGen.plotting.plot_wire_loops import plot_wire_loops, plot_wire_loops_tube + + +# ========================================================== +# Extract wire path with wire spacing +# ========================================================== +def extract_wire_paths( + verts, + faces, + psi, + levels, + display_debug_plots=False, + smooth_resample=4000, + wire_width=0.002, + wire_thickness=0.0016, + clearance=0.0005, + terminal_cut_length=0.0015): + + """ + Extract wire paths from a triangulated surface using contour levels. + This function generates spiral wire paths by extracting contours from a + triangulated surface, enforcing spacing constraints, and connecting loops + with smooth arc bridges. It handles both positive and negative polarity + coils separately. + Args: + verts (np.ndarray): Vertex coordinates of shape (n, 3) representing + the triangulated surface mesh. + faces (np.ndarray): Face connectivity array of shape (m, 3) defining + triangulation. + psi (np.ndarray): Scalar field values at vertices for contour extraction. + levels (array-like): Contour levels to extract from the scalar field. + display_debug_plots (bool, optional): If True, display visualization of + wire loops. Defaults to False. + smooth_resample (int, optional): Number of points for resampling smoothed + paths. Defaults to 4000. + wire_width (float, optional): Width of the wire in meters. Defaults to 0.002. + wire_thickness (float, optional): Thickness of the wire in meters. + Defaults to 0.0016. + clearance (float, optional): Minimum clearance between wires in meters. + Defaults to 0.0005. + terminal_cut_length (float, optional): Length to trim from terminal ends + in meters. Defaults to 0.0015. + Returns: + list: List of dictionaries, each containing: + - "path" (np.ndarray): 3D coordinates of wire centerline. + - "sign" (int): Coil polarity (1 for positive, -1 for negative). + - "width" (float): Wire width. + - "thickness" (float): Wire thickness. + Raises: + No explicit exceptions. Gracefully handles edge cases with fallbacks. + """ + z_plate = np.mean(verts[:, 2]) + + # Minimum spacing between wire centerlines + min_spacing = wire_width + clearance + 0.25 * wire_width + + triang = mtri.Triangulation(verts[:, 0], verts[:, 1], faces) + + fig, ax = plt.subplots() + cs = ax.tricontour(triang, psi, levels=levels) + plt.close(fig) + + # ------------------------------------------------ + # Smooth path + # ------------------------------------------------ + def smooth_path(path): + + if len(path) < 10: + return path + + try: + + k = min(3, len(path) - 1) + + tck, u = splprep( + [path[:,0], path[:,1], path[:,2]], + s=1e-6, + k=k + ) + + u_new = np.linspace(0,1,smooth_resample) + + x,y,z = splev(u_new,tck) + + return np.column_stack([x,y,z]) + + except: + return path + + # ------------------------------------------------ + # Arc bridge + # ------------------------------------------------ + def arc_bridge(p1, p2, z, npts=40): + + r1 = np.sqrt(p1[0]**2 + p1[1]**2) + r2 = np.sqrt(p2[0]**2 + p2[1]**2) + + r = 0.5 * (r1 + r2) + + th1 = np.arctan2(p1[1], p1[0]) + th2 = np.arctan2(p2[1], p2[0]) + + dth = th2 - th1 + + if dth > np.pi: + dth -= 2*np.pi + if dth < -np.pi: + dth += 2*np.pi + + th = np.linspace(th1, th1 + dth, npts) + + x = r * np.cos(th) + y = r * np.sin(th) + + return np.column_stack([x, y, np.full(npts, z)]) + + # ------------------------------------------------ + # Cut terminal end (NEW) + # ------------------------------------------------ + def cut_terminal_end(path, cut_length): + + if len(path) < 5: + return path + + seg_lengths = np.linalg.norm(np.diff(path, axis=0), axis=1) + + total = 0.0 + idx = len(path) - 1 + + while idx > 1 and total < cut_length: + total += seg_lengths[idx-1] + idx -= 1 + + return path[:idx] + + # ------------------------------------------------ + # Extract loops + # ------------------------------------------------ + loops = [] + + for level, segs in zip(cs.levels, cs.allsegs): + + for seg in segs: + + if len(seg) < 30: + continue + + r = np.mean(np.sqrt(seg[:,0]**2 + seg[:,1]**2)) + + path = np.column_stack([ + seg[:,0], + seg[:,1], + np.full(len(seg), z_plate) + ]) + + loops.append({ + "path": path, + "sign": 1 if level > 0 else -1, + "r": r + }) + + # Separate polarities + pos_loops = sorted([l for l in loops if l["sign"] > 0], key=lambda x: x["r"]) + neg_loops = sorted([l for l in loops if l["sign"] < 0], key=lambda x: x["r"]) + + # ------------------------------------------------ + # Enforce global spacing + # ------------------------------------------------ + def enforce_spacing(loopset): + + accepted = [] + global_pts = None + + for loop in loopset: + + pts = loop["path"] + + if global_pts is None: + + accepted.append(pts) + global_pts = pts + + continue + + tree = cKDTree(global_pts) + + d,_ = tree.query(pts) + + if np.min(d) > min_spacing: + + accepted.append(pts) + global_pts = np.vstack([global_pts, pts]) + + return accepted + + pos_loops = enforce_spacing(pos_loops) + neg_loops = enforce_spacing(neg_loops) + + # ------------------------------------------------ + # Build spiral + # ------------------------------------------------ + def build_spiral(loopset): + + if len(loopset) == 0: + return None + + spiral = loopset[0].copy() + + cut_pts = max(5, int(0.16 * len(spiral))) + spiral = spiral[:-cut_pts] + + for loop in loopset[1:]: + + prev_end = spiral[-1] + + tree = cKDTree(loop) + + _, idx = tree.query(prev_end) + + loop = np.roll(loop, -idx, axis=0) + + bridge = arc_bridge(prev_end, loop[0], z_plate) + + spiral = np.vstack([spiral, bridge, loop]) + + spiral = smooth_path(spiral) + + # NEW: trim terminal end + spiral = cut_terminal_end(spiral, terminal_cut_length) + + return spiral + + wire_loops = [] + + pos_spiral = build_spiral(pos_loops) + neg_spiral = build_spiral(neg_loops) + + if pos_spiral is not None: + wire_loops.append({ + "path": pos_spiral, + "sign": 1, + "width": wire_width, + "thickness": wire_thickness + }) + + if neg_spiral is not None: + wire_loops.append({ + "path": neg_spiral, + "sign": -1, + "width": wire_width, + "thickness": wire_thickness + }) + + # ------------------------------------------------ + # Visualization + # ------------------------------------------------ + if display_debug_plots: + plot_wire_loops_tube(wire_loops, wire_width) + + return wire_loops + diff --git a/pyCoilGen/sub_functions/generate_coil_from_stream_function.py b/pyCoilGen/sub_functions/generate_coil_from_stream_function.py new file mode 100644 index 00000000..3abf9993 --- /dev/null +++ b/pyCoilGen/sub_functions/generate_coil_from_stream_function.py @@ -0,0 +1,393 @@ +import numpy as np +import trimesh +import matplotlib.pyplot as plt +from scipy.spatial import cKDTree +from scipy.ndimage import gaussian_filter1d +from shapely.geometry import Point +from collections import defaultdict +from mpl_toolkits.mplot3d import Axes3D +############################################################ +# PARAMETERS +############################################################ + +N_LEVELS = 12 +WIRE_RADIUS = 0.001 # 1 mm channel radius (~AWG14) +SMOOTH_SIGMA = 1.0 + +############################################################ +# 1 EXTRACT STREAM FUNCTION +############################################################ + +def extract_stream_function(coil_parts): + + meshes = [] + psis = [] + + for part in coil_parts: + meshes.append(part.coil_mesh) + psis.append(part.stream_function) + + + + return meshes, psis + + +############################################################ +# 2 NORMALIZE STREAM FUNCTION +############################################################ + +def normalize_stream_function(psi): + + psi = psi - np.mean(psi) + psi = psi / np.max(np.abs(psi)) + + return psi + + +############################################################ +# 3 ROBUST CONTOUR LEVELS +############################################################ + +def compute_levels(psi, n_levels=N_LEVELS): + + levels = np.quantile( + psi, + np.linspace(0.05, 0.95, n_levels) + ) + + return levels + + +############################################################ +# 4 TRIANGLE CONTOUR INTERSECTION +############################################################ + +def triangle_contour(vertices, psi_vals, level): + + edges = [(0,1),(1,2),(2,0)] + points = [] + + for i,j in edges: + + p1,p2 = psi_vals[i], psi_vals[j] + + if (p1-level)*(p2-level) < 0: + + t = (level-p1)/(p2-p1) + + point = vertices[i] + t*(vertices[j]-vertices[i]) + + points.append(point) + + if len(points) == 2: + return points + + return None + + +############################################################ +# 5 COMPUTE ALL CONTOUR SEGMENTS +############################################################ + + + +# --------------------------------------------------------- +# Helper: compute intersection of edge with contour level +# --------------------------------------------------------- + +def edge_intersection(v1, v2, p1, p2, level): + + if (p1 - level) * (p2 - level) > 0: + return None + + if abs(p1 - p2) < 1e-12: + return None + + t = (level - p1) / (p2 - p1) + + if t < 0 or t > 1: + return None + + return v1 + t * (v2 - v1) + + +# --------------------------------------------------------- +# Marching triangles +# --------------------------------------------------------- + +def triangle_contour(vertices, psi_vals, level): + + v0, v1, v2 = vertices + p0, p1, p2 = psi_vals + + pts = [] + + e = edge_intersection(v0, v1, p0, p1, level) + if e is not None: + pts.append(e) + + e = edge_intersection(v1, v2, p1, p2, level) + if e is not None: + pts.append(e) + + e = edge_intersection(v2, v0, p2, p0, level) + if e is not None: + pts.append(e) + + if len(pts) == 2: + return np.array(pts) + + return None + + +# --------------------------------------------------------- +# Segment stitching +# --------------------------------------------------------- + +def stitch_segments(segments, tol=20e-3): + + segments = [tuple(map(tuple, s)) for s in segments] + + unused = set(range(len(segments))) + loops = [] + + while unused: + + idx = unused.pop() + seg = segments[idx] + + loop = [np.array(seg[0]), np.array(seg[1])] + + growing = True + + while growing: + growing = False + + for j in list(unused): + + s = segments[j] + + p1 = np.array(s[0]) + p2 = np.array(s[1]) + + if np.linalg.norm(loop[-1] - p1) < tol: + loop.append(p2) + unused.remove(j) + growing = True + break + + if np.linalg.norm(loop[-1] - p2) < tol: + loop.append(p1) + unused.remove(j) + growing = True + break + + if np.linalg.norm(loop[0] - p1) < tol: + loop.insert(0, p2) + unused.remove(j) + growing = True + break + + if np.linalg.norm(loop[0] - p2) < tol: + loop.insert(0, p1) + unused.remove(j) + growing = True + break + + loops.append(np.array(loop)) + print(f"Formed loop with {len(loop)} points, remaining segments: {len(unused)}") + + return loops + + +# --------------------------------------------------------- +# Main contour extraction +# --------------------------------------------------------- + +def compute_streamfunction_loops(mesh, psi, levels, display_loops:bool = True): + + vertices = mesh.v + faces = mesh.f + + segments_by_level = defaultdict(list) + + for face in faces: + + tri_v = vertices[face] + tri_p = psi[face] + + for level in levels: + + seg = triangle_contour(tri_v, tri_p, level) + + if seg is not None: + segments_by_level[level].append(seg) + + all_loops = [] + + for level, segs in segments_by_level.items(): + + print("Level", level, "segments:", len(segs)) + + loops = stitch_segments(segs) + + print("Level", level, "loops:", len(loops)) + + # remove tiny fragments + loops = [l for l in loops if len(l) > 10] + + all_loops.extend(loops) + + print("Total loops:", len(all_loops)) + + if display_loops: + + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + + for loop in loops: + ax.plot(loop[:,0], loop[:,1], loop[:,2], 'r') + + ax.set_title("Extracted Gradient Coil Loops") + + plt.show() + + return all_loops + + + + +############################################################ +# 6 SMOOTH LOOPS +############################################################ + +def smooth_loops(loops): + + smoothed = [] + + for loop in loops: + + if len(loop) < 5: + smoothed.append(loop) + continue + + sm = gaussian_filter1d(loop, SMOOTH_SIGMA, axis=0) + + smoothed.append(sm) + + return smoothed + + +############################################################ +# 8 VISUALIZE LOOPS +############################################################ + +def visualize_loops(loops): + + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + + for loop in loops: + + ax.plot(loop[:,0], loop[:,1], loop[:,2]) + + ax.set_title("Gradient Coil Wire Patterns") + + plt.show() + + +############################################################ +# 9 CREATE PRINTABLE WIRE CHANNELS +############################################################ + +def create_wire_channels(loops, radius=WIRE_RADIUS): + + tubes = [] + profile = Point(0, 0).buffer(radius) + + for loop in loops: + + tube = trimesh.creation.sweep_polygon( + polygon=profile, + path=loop, + cap=False # Optional: caps the ends + ) + + tubes.append(tube) + + return trimesh.util.concatenate(tubes) + +############################################################# +# Visualize the stream function values on the mesh surface +############################################################# + +def plot_stream_function(meshes, psis, cmap='coolwarm'): + fig = plt.figure(figsize=(12, 6)) + + for i, (mesh, psi) in enumerate(zip(meshes, psis)): + ax = fig.add_subplot(1, len(meshes), i+1, projection='3d') + + # Get vertices and faces + verts = mesh.v + faces = mesh.f + + # Create a triangulated surface + tri_collection = ax.plot_trisurf(verts[:,0], verts[:,1], verts[:,2], + triangles=faces, + cmap=cmap, + linewidth=0.2, + antialiased=True, + shade=True) + + # Map stream function to color + tri_collection.set_array(psi) + tri_collection.autoscale() + + ax.set_title(f'Coil part {i+1}') + ax.set_axis_off() + + fig.colorbar(tri_collection, ax=fig.axes, shrink=0.5, label='Stream Function') + plt.show() + + + +############################################################ +# 10 COMPLETE PIPELINE +############################################################ + +def generate_coil_from_stream_function(coil_parts, n_levels, display=True): + + meshes, psis = extract_stream_function(coil_parts) + + if display: + plot_stream_function(meshes, psis) + + all_loops = [] + + for mesh, psi in zip(meshes, psis): + + print("Processing mesh with", len(mesh.v), "vertices") + + psi = normalize_stream_function(psi) + + levels = compute_levels(psi, n_levels) + + loops = compute_streamfunction_loops(mesh, psi, levels) + + loops = smooth_loops(loops) + + all_loops.extend(loops) + + print("Total loops:", len(all_loops)) + + visualize_loops(all_loops) + + channels = create_wire_channels(all_loops) + + channels.export("gradient_wire_channels.stl") + + print("Exported gradient_wire_channels.stl") + + return all_loops, channels + + diff --git a/pyCoilGen/sub_functions/get_streamfunction.py b/pyCoilGen/sub_functions/get_streamfunction.py new file mode 100644 index 00000000..95815770 --- /dev/null +++ b/pyCoilGen/sub_functions/get_streamfunction.py @@ -0,0 +1,133 @@ +import numpy as np +from scipy.sparse.linalg import lsmr + + +def solve_streamfunction_with_initial_guess( + reduced_sensitivity_matrix, + target_field_single, + mesh_vertices, + gradient_axis="x", + prior_strength=0.0, + atol=1e-9, + btol=1e-9, + maxiter=1000): + + """ + Solve for streamfunction coefficients using LSMR with a physics-informed prior. + + Parameters + ---------- + reduced_sensitivity_matrix : ndarray (M x K) + Sensitivity matrix in reduced basis. + + target_field_single : ndarray (M,) + Target field values. + + mesh_vertices : ndarray (N x 3) + Coil mesh vertices. + + streamfunction_basis : ndarray (N x K) + Basis matrix mapping coefficients -> streamfunction. + + gradient_axis : str + 'x', 'y', or 'z'. + + prior_strength : float + Weight of the prior regularization. + + Returns + ------- + psi : ndarray (N,) + Streamfunction on mesh vertices. + + coeffs : ndarray (K,) + Reduced basis coefficients. + + B_pred : ndarray (M,) + Predicted magnetic field. + """ + + A = reduced_sensitivity_matrix + b = target_field_single + verts = mesh_vertices + n_verts = verts.shape[0] + streamfunction_basis = np.eye(n_verts) + + # ----------------------------------------- + # 1. Build physics-informed ψ prior + # ----------------------------------------- + + x = verts[:,0] + y = verts[:,1] + z = verts[:,2] + + + + if gradient_axis == "x": + psi0 = y + elif gradient_axis == "y": + psi0 = -x + elif gradient_axis == "z": + psi0 = x * y + else: + raise ValueError("gradient_axis must be x, y, or z") + + # normalize prior + psi0 = psi0 - np.mean(psi0) + psi0 = psi0 / np.max(np.abs(psi0)) + + # ----------------------------------------- + # 2. Project ψ prior into reduced basis + # ----------------------------------------- + + x0, *_ = np.linalg.lstsq(streamfunction_basis, psi0, rcond=None) + + # ----------------------------------------- + # 3. Apply prior regularization + # ----------------------------------------- + + if prior_strength > 0: + + A_aug = np.vstack([ + A, + prior_strength * np.eye(A.shape[1]) + ]) + + b_aug = np.concatenate([ + b, + prior_strength * x0 + ]) + + else: + + A_aug = A + b_aug = b + + # ----------------------------------------- + # 4. Solve with LSMR + # ----------------------------------------- + + result = lsmr( + A_aug, + b_aug, + x0=x0, + atol=atol, + btol=btol, + maxiter=maxiter + ) + + coeffs = result[0] + + # ----------------------------------------- + # 5. Reconstruct streamfunction + # ----------------------------------------- + + psi = streamfunction_basis @ coeffs + + # ----------------------------------------- + # 6. Compute predicted field + # ----------------------------------------- + + B_pred = A @ coeffs + + return psi, coeffs, B_pred \ No newline at end of file diff --git a/pyCoilGen/sub_functions/gradient_former.py b/pyCoilGen/sub_functions/gradient_former.py new file mode 100644 index 00000000..13b7cb06 --- /dev/null +++ b/pyCoilGen/sub_functions/gradient_former.py @@ -0,0 +1,409 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle, Circle, FancyArrow +import trimesh +from scipy.interpolate import splprep, splev +from pyCoilGen.plotting.plot_wire_loops import plot_gerber_paths, plot_stl_patch + +# ========================================================== +# 3D printer G-code generation for gradient coil holder +# ========================================================== + +def generate_gradient_former( + + coil_parts, + output_prefix="gradient_plate", + plate_diameter_m=0.1524, + plate_thickness_m=0.003, + wire_width_m=0.002, + display_debug_plots=True, + save_loop_coords=True): + + """ + Generate gradient former coil plates with Gerber files and STL models. + This function processes coil parts to generate PCB gradient coils with + traces, plate geometries, and 3D printable STL files. It creates Gerber + files for manufacturing and optional 3D models with features including + tubes, alignment pins, terminal holes, and strain relief. + Args: + coil_parts (list): List of coil part objects, each containing wire_loops + with path and sign attributes. + output_prefix (str, optional): Prefix for output filenames. + Defaults to "gradient_plate". + plate_diameter_m (float, optional): Diameter of the plate in meters. + Defaults to 0.1524 (6 inches). + plate_thickness_m (float, optional): Thickness of the plate in meters. + Defaults to 0.003 (3 mm). + wire_width_m (float, optional): Width of wire traces in meters. + Defaults to 0.002 (2 mm). + display_debug_plots (bool, optional): Whether to display debug + visualization plots. Defaults to True. + save_loop_coords (bool, optional): Whether to save loop coordinates + to text files. Defaults to True. + Returns: + None + Raises: + Exception: If STL generation or boolean operations fail during + 3D model creation. + Note: + Generated files include: + - Gerber track files (.gbr) + - Complement Gerber files (.gbr) + - Loop coordinate files (.txt) if save_loop_coords is True + - Final STL models (.stl) for 3D printing + """ + + + + + mm = 1000.0 + plate_radius_mm = plate_diameter_m*mm/2 + wire_width_mm = wire_width_m*mm + + #------------------------- + # Generate the plate STL for merging later if needed + #------------------------- + + + + # ---------------------------------------- + # Helper to ensure minimum spacing between points + # ---------------------------------------- + def enforce_global_spacing(points, min_dist_mm=0.1): + filtered=[points[0]] + for p in points[1:]: + if np.linalg.norm(p-filtered[-1])>min_dist_mm: + filtered.append(p) + return np.array(filtered) + + # ---------------------------------------- + # Helper to generate circular plate + # ---------------------------------------- + def generate_plate_circle(radius_mm, npts=720): + th = np.linspace(0,2*np.pi,npts) + x = radius_mm*np.cos(th) + y = radius_mm*np.sin(th) + return np.column_stack([x,y]) + + # ---------------------------------------- + # Process each coil plate + # ---------------------------------------- + for part_index, part in enumerate(coil_parts): + + print(f"\nProcessing plate {part_index}") + + gerber_paths = [] + + # ---------------------------------------- + # Collect paths + # ---------------------------------------- + for loop_index, loop in enumerate(part.wire_loops): + + pts_mm = np.asarray(loop["path"])[:, :2]*mm + pts_mm = enforce_global_spacing(pts_mm, min_dist_mm=0.5) + + gerber_paths.append((pts_mm, loop["sign"])) + + if save_loop_coords: + coords = np.column_stack([ + pts_mm[:,0], + pts_mm[:,1], + np.zeros(len(pts_mm)) + ]) + + fname = f"{output_prefix}_plate{part_index}_loop{loop_index}.txt" + np.savetxt(fname, coords, fmt="%.2f") + + + + # ---------------------------------------- + # Write standard Gerber (tracks) + # ---------------------------------------- + gerber_fname = f"{output_prefix}_plate{part_index}.gbr" + + with open(gerber_fname, "w") as f: + + f.write("G04 Gradient coil generated by Python*\n") + f.write("%MOMM*%\n") + f.write("%FSLAX24Y24*%\n") + + f.write(f"%ADD10C,{wire_width_mm:.4f}*%\n") + + f.write("D10*\n") + + for path, sign in gerber_paths: + + x0, y0 = path[0] + + f.write(f"X{int(x0*10000):07d}Y{int(y0*10000):07d}D02*\n") + + for p in path[1:]: + x, y = p + f.write(f"X{int(x*10000):07d}Y{int(y*10000):07d}D01*\n") + + f.write("M02*\n") + + print("Saved Gerber:", gerber_fname) + + # ---------------------------------------- + # Write complement Gerber (plate minus tracks) + # ---------------------------------------- + complement_fname = f"{output_prefix}_plate{part_index}_complement.gbr" + + plate_circle = generate_plate_circle(plate_radius_mm) + + with open(complement_fname, "w") as f: + + f.write("G04 Complement copper layer*\n") + f.write("%MOMM*%\n") + f.write("%FSLAX24Y24*%\n") + + # aperture for tracks + f.write(f"%ADD10C,{wire_width_mm:.4f}*%\n") + + # ---- Filled plate region ---- + f.write("%LPD*%\n") + f.write("G36*\n") + + x0, y0 = plate_circle[0] + + f.write(f"X{int(x0*10000):07d}Y{int(y0*10000):07d}D02*\n") + + for p in plate_circle[1:]: + x,y = p + f.write(f"X{int(x*10000):07d}Y{int(y*10000):07d}D01*\n") + + f.write(f"X{int(x0*10000):07d}Y{int(y0*10000):07d}D01*\n") + + f.write("G37*\n") + + # ---- Subtract wire paths ---- + f.write("%LPC*%\n") + f.write("D10*\n") + + for path,_ in gerber_paths: + + x0,y0 = path[0] + + f.write(f"X{int(x0*10000):07d}Y{int(y0*10000):07d}D02*\n") + + for p in path[1:]: + x,y = p + f.write(f"X{int(x*10000):07d}Y{int(y*10000):07d}D01*\n") + + f.write("M02*\n") + + print("Saved Complement Gerber:", complement_fname) + + # ---------------------------------------- + # Visualization + # ---------------------------------------- + if display_debug_plots: + plot_gerber_paths(gerber_paths, plate_radius_mm, wire_width_mm, part_index) + # ---------------------------------------- + # STL GENERATION (FINAL: all features) + # ---------------------------------------- + try: + + + def smooth_path(x, y, smoothing=5.0, n_points=300): + tck, _ = splprep([x, y], s=smoothing, per=True) + u = np.linspace(0, 1, n_points) + x_new, y_new = splev(u, tck) + return np.array(x_new), np.array(y_new) + + all_segments = [] + hole_meshes = [] + strain_relief_segments = [] + numbering_meshes = [] + + plate_height_mm = plate_thickness_m * mm + hole_radius = wire_width_mm * 0.8 + hole_height = plate_height_mm * 2.0 + offset_dist = wire_width_mm * 1.5 + strain_length = wire_width_mm * 4.0 + tube_radius = wire_width_mm / 2.0 + + # ---- Plate ---- + plate = trimesh.creation.cylinder( + radius=plate_radius_mm, + height=plate_height_mm, + sections=128 + ) + plate.apply_translation([0, 0, plate_height_mm / 2]) + + # ---- Alignment pins (snap-fit) ---- + pin_radius = wire_width_mm + pin_height = plate_height_mm * 1.2 + pin_offset = plate_radius_mm * 0.85 + pin_angles = [0, np.pi/2, np.pi, 3*np.pi/2] # 4 pins + + pin_meshes = [] + for a in pin_angles: + x = pin_offset * np.cos(a) + y = pin_offset * np.sin(a) + pin = trimesh.creation.cylinder( + radius=pin_radius, + height=pin_height + ) + pin.apply_translation([x, y, plate_height_mm / 2]) + pin_meshes.append(pin) + + # ---- Process each loop ---- + for loop_index, (path, sign) in enumerate(gerber_paths): + + x_raw, y_raw = path[:,0], path[:,1] + points = smooth_path(x_raw, y_raw, smoothing=5.0) + + points_3d = np.column_stack([points[0], points[1], np.zeros_like(points[0])]) + + # ---- Tube segments ---- + for i in range(len(points_3d)-1): + p1, p2 = points_3d[i], points_3d[i+1] + if np.linalg.norm(p2 - p1) < 0.5: + continue + seg = trimesh.creation.cylinder( + radius=tube_radius, + segment=[p1, p2] + ) + # optional multi-layer: raise every other loop slightly + if loop_index % 2 == 1: + seg.apply_translation([0,0,0.2]) + all_segments.append(seg) + + # ---- Terminal holes ---- + def tangent(p0, p1): + v = p1 - p0 + n = np.linalg.norm(v) + return v/n if n>0 else np.array([1.0,0.0]) + def perp(v): return np.array([-v[1], v[0]]) + + t_start = tangent(path[0], path[1]) + t_end = tangent(path[-2], path[-1]) + start_pt = path[0] + perp(t_start)*offset_dist + end_pt = path[-1] + perp(t_end)*offset_dist + + for pt in [start_pt, end_pt]: + hole = trimesh.creation.cylinder( + radius=hole_radius, + height=hole_height + ) + hole.apply_translation([pt[0], pt[1], plate_height_mm / 2]) + hole_meshes.append(hole) + + # ---- Strain relief ---- + for base_pt, t_vec in [(start_pt, t_start), (end_pt, t_end)]: + p1 = np.array([base_pt[0], base_pt[1], 0]) + p2 = p1 + np.append(t_vec * strain_length, 0) + relief = trimesh.creation.cylinder( + radius=tube_radius, + segment=[p1, p2] + ) + strain_relief_segments.append(relief) + + # ---- Loop numbering ---- + # try: + # num_text = trimesh.creation.text( + # text=str(loop_index+1), + # height=wire_width_mm*2, + # depth=0.5 + # ) + # # place number at first point of loop, slightly above plate + # num_text.apply_translation([ + # path[0][0], + # path[0][1], + # plate_height_mm - 0.3 + # ]) + # numbering_meshes.append(num_text) + # except Exception as e: + # print(f"⚠️ Loop numbering failed for loop {loop_index+1}: {e}") + + # ---- Combine all cutters ---- + tube = trimesh.util.concatenate(all_segments) + if strain_relief_segments: + relief_mesh = trimesh.util.concatenate(strain_relief_segments) + tube = trimesh.util.concatenate([tube, relief_mesh]) + cutter_meshes = [tube] + hole_meshes + numbering_meshes + pin_meshes + cutter = trimesh.util.concatenate(cutter_meshes) + + # ---- Offset tube slightly into plate for engraving/groove ---- + cutter.apply_translation([0,0,plate_height_mm*0.75]) + + # # ---- Plate engraving ---- + # label = "TOP" if part_index==0 else "BOTTOM" + # try: + # text_mesh = trimesh.creation.text( + # text=label, + # height=plate_radius_mm*0.15, + # depth=0.5 + # ) + # text_mesh.apply_translation([ + # -plate_radius_mm*0.3, + # -plate_radius_mm*0.8, + # plate_height_mm - 0.5 + # ]) + # cutter = trimesh.util.concatenate([cutter, text_mesh]) + # except: + # print("⚠️ Plate engraving skipped") + + # ---- Boolean difference ---- + result = plate.difference(cutter) + if result is None or len(result.vertices)==0: + print("❌ Boolean failed") + continue + + # ---- Export STL ---- + stl_fname = f"{output_prefix}_plate{part_index}_final_print.stl" + result.export(stl_fname) + print(f"✅ Final STL saved: {stl_fname}") + + + # ---------------------------------------- + # DEBUG VISUALIZATION (Matplotlib + terminals + strain relief + loop numbers) + # ---------------------------------------- + if display_debug_plots: + plot_stl_patch(gerber_paths, plate_radius_mm, wire_width_mm, part_index) + + except Exception as e: + print("❌ STL generation failed:", e) + + print("\nGradient former coordinate export complete with Gerber and visualization.") + + + + +# =============================== +# CHECK SPHERICAL HARMONICS +#================================ + + +def run_spherical_harmonic_diagnostic(result): + + coords = result.target_field.coords + field = result.coil_gradient.gradient_in_target_direction + + x,y,z = coords + + r = np.sqrt(x**2+y**2+z**2) + theta = np.arccos(z/r) + phi = np.arctan2(y,x) + + max_l = 3 + + coeffs = {} + + for l in range(max_l+1): + for m in range(-l,l+1): + + Ylm = sph_harm_y(l,m,phi,theta).real + a = np.sum(field*Ylm)/np.sum(Ylm**2) + + coeffs[(l,m)] = a + + print("\n===== Spherical Harmonic Content =====") + + for k,v in coeffs.items(): + print(f"l={k[0]}, m={k[1]} : {v:.4e}") + + + diff --git a/pyCoilGen/sub_functions/parse_input.py b/pyCoilGen/sub_functions/parse_input.py index bc326e75..b0165a43 100644 --- a/pyCoilGen/sub_functions/parse_input.py +++ b/pyCoilGen/sub_functions/parse_input.py @@ -189,6 +189,9 @@ def parse_input(parse_cli=True): parser.add_argument('--smooth_factor', type=int, default=1, help="Smoothing parameter if tracks should be smoothed (i.e. when > 1).") + # Add flag for cylindrical surface + parser.add_argument('--target_field_weighting', type=bool, default=False, + help="Increased weighting in the center of the target field region, which can lead to improved gradient efficiency at the cost of homogeneity") """ Unused # Add flag to plot results parser.add_argument('--plot_flag', type=bool, default=True, @@ -317,6 +320,8 @@ def parse_input(parse_cli=True): # Add the parameters for the generation of the (default) biplanar mesh parser.add_argument('--persistence_dir', type=str, default='debug', help=f"Directory to write persistence data") + + # Add a parameter for the location of the FastHenry2 executable, platform aware os_system = platform.system() if os_system == "Windows": diff --git a/pyCoilGen/sub_functions/simulate_gradient_coil.py b/pyCoilGen/sub_functions/simulate_gradient_coil.py new file mode 100644 index 00000000..44742ee9 --- /dev/null +++ b/pyCoilGen/sub_functions/simulate_gradient_coil.py @@ -0,0 +1,91 @@ +import numpy as np +import magpylib as magpy +import matplotlib.pyplot as plt +from pyCoilGen.plotting.plot_field import plot_field +import pyCoilGen.plotting.plot_wire_loops as plot_wire_loops + +# ========================================================== +# Simulate using magpylib +# ========================================================== +def simulate_gradient_coil(coil_parts, DSV_IMAGING, wire, display_field=False): + """ + Simulate gradient coil B-field including wire width and thickness. + + Parameters + ---------- + coil_parts : list + List of coil parts, each having `wire_loops` extracted from extract_wire_paths. + DSV_IMAGING : float + Diameter of the imaging sphere in meters. + wire : dict + Contains {"current": } in Amperes. + display_field : bool + If True, display 3D scatter plots of Bx, By, Bz. + + Returns + ------- + dict + {"points": points, "B": B, "Bx": Bx, "By": By, "Bz": Bz} + """ + + + # ----------------------------- + # Create imaging sphere points + # ----------------------------- + r = DSV_IMAGING / 2 + grid = np.linspace(-r, r, 21) + points = np.array([[x, y, z] + for x in grid for y in grid for z in grid + if x**2 + y**2 + z**2 <= r**2]) + + current_scale = wire["current"] + B = np.zeros((len(points), 3)) + + # ----------------------------- + # Loop over coil parts and loops + # ----------------------------- + for part in coil_parts: + for spiral in part.wire_loops: + path = spiral["path"] + width = spiral.get("width", wire['width']) # default to thickness if width not specified + thickness = spiral.get("thickness", wire['thickness']) # default to width if thickness not specified + + # Downsample if spline too long (memory efficiency) + if len(path) > 400: + idx = np.linspace(0, len(path)-1, 400).astype(int) + path = path[idx] + + # 4-point rectangular quadrature to approximate wire cross-section + offsets = [(-width/2, -thickness/2), + (-width/2, thickness/2), + ( width/2, -thickness/2), + ( width/2, thickness/2)] + + for dx, dz in offsets: + offset_path = path.copy() + offset_path[:, 0] += dx + offset_path[:, 2] += dz + + # Polyline source (current direction taken from vertex order) + source = magpy.current.Polyline( + current=current_scale / len(offsets), # divided among offsets + vertices=offset_path + ) + + # accumulate B-field + B += magpy.getB(source, points) + + # ----------------------------- + # Extract components + # ----------------------------- + Bx = B[:, 0] + By = B[:, 1] + Bz = B[:, 2] + + # ----------------------------- + # Optional visualization + # ----------------------------- + if display_field: + plot_field(Bx, By, Bz, points) + + return {"points": points, "B": B, "Bx": Bx, "By": By, "Bz": Bz} diff --git a/pyCoilGen/sub_functions/stl_mesh_generation.py b/pyCoilGen/sub_functions/stl_mesh_generation.py new file mode 100644 index 00000000..d11956b4 --- /dev/null +++ b/pyCoilGen/sub_functions/stl_mesh_generation.py @@ -0,0 +1,217 @@ +import numpy as np +import struct +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + +radius = 0.07 +nr = 30 # number of radial divisions +nt = 90 # 180 for 2 degree resolution, 90 for 4 degree resolution +inner_radius = radius / nr +z_positions = [-0.03675, 0.03675] +filename = "data/pyCoilGenData/Geometry_Data/Tenacity_circular_1.stl" + + +def create_mesh(radius,nr,nt,z): + """ + Creates a mesh for a cylindrical surface. + Generates vertices and faces for a mesh representing a cylindrical surface + with a specified radius and height. The mesh is created by arranging vertices + in concentric circles and connecting them with triangular faces. + Args: + radius (float): The outer radius of the cylindrical mesh. + nr (int): The number of radial divisions (rings). + nt (int): The number of theta (angular) divisions (sectors). + z (float): The z-coordinate (height) of the mesh plane. + Returns: + tuple: A tuple containing: + - verts (np.ndarray): An (nr*nt, 3) array of vertex coordinates [x, y, z]. + - faces (np.ndarray): An (n, 3) array of triangular face indices. + Note: + This function assumes `inner_radius` is defined in the outer scope. + The mesh is created as a flat cylindrical annulus at height z. + """ + + + verts=[] + faces=[] + + for i in range(nr): + + r = inner_radius + (radius-inner_radius)*i/(nr-1) + + for j in range(nt): + + theta=2*np.pi*j/nt + + x=r*np.cos(theta) + y=r*np.sin(theta) + + verts.append([x,y,z]) + + verts=np.array(verts) + + def vid(i,j): + return i*nt+j + + for i in range(nr-1): + + for j in range(nt): + + v1=vid(i,j) + v2=vid(i,(j+1)%nt) + v3=vid(i+1,(j+1)%nt) + v4=vid(i+1,j) + + faces.append([v1,v2,v3]) + faces.append([v1,v3,v4]) + + return verts,np.array(faces) + + +def write_stl(filename,verts,faces): + """ + Write a mesh to an STL (Stereolithography) file format. + This function takes a set of vertices and faces defining a 3D mesh and writes + them to a binary STL file. It ensures that all triangle normals point outward + by flipping triangles as needed based on the Z-component of the normal vector. + Args: + filename (str): The path and name of the output STL file to be created. + verts (numpy.ndarray): An array of vertices with shape (N, 3) containing + the 3D coordinates of each vertex in the mesh. + faces (list or numpy.ndarray): A list or array of triangular faces where + each face is defined by three indices referring to vertices in the + verts array. + Returns: + None: The function writes directly to a file instead of returning a value. + Notes: + - The function modifies the faces array in place to ensure consistent + triangle winding order. + - Triangle normals with negative Z-components are flipped to ensure + outward-facing normals. + - Degenerate triangles (with zero area) have their normal set to [0, 0, 0]. + - The output file uses binary STL format with an 80-byte header and + 4-byte face count. + Raises: + IOError: If the file cannot be opened or written to. + ValueError: If faces contain invalid vertex indices. + """ + + + for i,tri in enumerate(faces): + + v1,v2,v3 = verts[tri] + + normal = np.cross(v2-v1, v3-v1) + + if normal[2] < 0: # flip triangle + faces[i] = [tri[0], tri[2], tri[1]] + + with open(filename,"wb") as f: + + f.write(b'PythonSTL'+b' '*(80-len('PythonSTL'))) + f.write(struct.pack(">> create_stl_mesh(radius=10, nr=20, nt=32, + ... z_positions=[0, 5, 10], filename='mesh.stl') + STL saved: mesh.stl + """ + + all_verts=[] + all_faces=[] + offset=0 + for z in z_positions: + + v,f=create_mesh(radius,nr,nt,z) + + all_verts.append(v) + all_faces.append(f+offset) + + offset+=len(v) + + verts=np.vstack(all_verts) + faces=np.vstack(all_faces) + + write_stl(filename,verts,faces) + + print("STL saved:",filename) \ No newline at end of file diff --git a/pyCoilGen/sub_functions/stream_function_optimization.py b/pyCoilGen/sub_functions/stream_function_optimization.py index 8b50194d..0f2bfb30 100644 --- a/pyCoilGen/sub_functions/stream_function_optimization.py +++ b/pyCoilGen/sub_functions/stream_function_optimization.py @@ -1,7 +1,7 @@ import numpy as np from scipy.optimize import minimize import ast # For minimization function parameter parsing - +from scipy.sparse.linalg import lsqr from typing import List # Logging @@ -11,11 +11,16 @@ from .data_structures import DataStructure, CoilPart from .constants import get_level, DEBUG_NONE, DEBUG_VERBOSE from pyCoilGen.helpers.common import blkdiag +# from .get_streamfunction import solve_streamfunction_with_initial_guess +from scipy.sparse import csr_matrix, diags +from scipy.sparse.linalg import lsqr +from matplotlib.tri import Triangulation +# from scipy.sparse import csr_matrix, block_diag log = logging.getLogger(__name__) -def stream_function_optimization(coil_parts: List[CoilPart], target_field, input_args) -> (List[CoilPart], np.ndarray, np.ndarray): +def stream_function_optimization(coil_parts: List[CoilPart], target_field, input_args) -> tuple[List[CoilPart], np.ndarray, np.ndarray]: """ Performs stream function optimization on coil parts. @@ -140,13 +145,11 @@ def stream_function_optimization(coil_parts: List[CoilPart], target_field, input else: # For initialization, calculate the Tikhonov solution; then do an iterative optimization log.info("Optimising with %s function.", input_args.sf_opt_method) - tik_reg_mat = tikhonov_reg_factor * reduced_res_matrix - reduced_sf = np.linalg.pinv( - reduced_sensitivity_matrix.T - @ reduced_sensitivity_matrix - + tik_reg_mat.T - @ tik_reg_mat - ) @ reduced_sensitivity_matrix.T @ target_field_single + + # Initial estimate using lsqr instead of pinv - + A = reduced_sensitivity_matrix + b = target_field_single + reduced_sf = lsqr(A, b, atol=1e-9, btol=1e-9, iter_lim=1000)[0] # Find the constrained solution stream_func_max = np.max(reduced_sf) * 2 @@ -154,9 +157,27 @@ def stream_function_optimization(coil_parts: List[CoilPart], target_field, input ub = np.ones_like(reduced_sf) * stream_func_max def cost_function(x): - return np.sum( - (reduced_sensitivity_matrix @ x - target_field_single) ** 2 - ) + tikhonov_reg_factor * (x.T @ reduced_res_matrix @ x) + # beta0 = 1e3 + # return np.sum( + # ((reduced_sensitivity_matrix @ x - target_field_single) ** 2 + # ) + tikhonov_reg_factor * (x.T @ reduced_res_matrix @ x) + beta0 * np.mean(reduced_sensitivity_matrix @ x)** 2) + alpha0 = 1e2 + beta0 = 1e2 + beta_z = 1e2 + B = reduced_sensitivity_matrix @ x + mean_field = np.mean(B) + data_term = np.sum((B - target_field_single) ** 2) + reg_term = x.T @ reduced_res_matrix @ x + # z_mode = z_coords / np.max(np.abs(z_coords)) + # coeff_z = np.dot(B, z_mode) / np.dot(z_mode, z_mode) + + cost = alpha0 * data_term \ + + tikhonov_reg_factor * (reg_term) \ + + beta0 * mean_field**2 + + print(data_term, reg_term, mean_field) + return cost + # Define the bounds bounds = [(lb[i], ub[i]) for i in range(len(lb))] @@ -179,18 +200,37 @@ def cost_function(x): input_args.minimize_method, method_params, minimize_method_options) # Perform the optimization - result = minimize(cost_function, reduced_sf, method=input_args.minimize_method, bounds=bounds, - **method_params, options=minimize_method_options) - log.debug("Detailed minimize result: Success: %s, message: %s", result.success, result.message) - if get_level() >= DEBUG_VERBOSE: - log.debug("Detailed minimize result: %s", result) - reduced_sf = result.x + def callback(xk): + if callback.count % 1 == 0: + log.info("Iteration %d, Cost: %.6e", callback.count, cost_function(xk)) + callback.count += 1 + + callback.count = 0 + # result = minimize(cost_function, reduced_sf, method=input_args.minimize_method, bounds=bounds, + # callback=callback, **method_params, options=minimize_method_options) + # log.debug("Detailed minimize result: Success: %s, message: %s", result.success, result.message) + + A = reduced_sensitivity_matrix + R = reduced_res_matrix + b = target_field_single + + ATA = A.T @ A + ATb = A.T @ b + + reduced_sf = np.linalg.solve(ATA + tikhonov_reg_factor * R, ATb) + + + + # if get_level() >= DEBUG_VERBOSE: + # log.debug("Detailed minimize result: %s", result) + # reduced_sf = result.x # Reexpand the stream potential to the boundary nodes opt_stream_func = reexpand_stream_function_for_boundary_nodes( reduced_sf, boundary_nodes, is_not_boundary_node, set_zero_flag) combined_mesh.stream_function = opt_stream_func + # Calculate the magnetic field generated by the optimized stream function b_field_opt_sf = np.vstack( ( @@ -200,6 +240,24 @@ def cost_function(x): ) ).T + # Visualize the magnetic field generated by the optimized stream function in the target region + if get_level() >= DEBUG_VERBOSE: + import matplotlib.pyplot as plt + coords = target_field.coords + x = coords[0, :] + y = coords[1, :] + z = coords[2, :] + b_field_opt_sf_x = b_field_opt_sf[:, 0] + b_field_opt_sf_y = b_field_opt_sf[:, 1] + b_field_opt_sf_z = b_field_opt_sf[:, 2] + + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + sc = ax.scatter(x, y, z, c=b_field_opt_sf_z * 1e3, cmap="jet") + plt.colorbar(sc, label="B-field (mT)") + ax.set_title("B-field generated by optimized stream function in target region") + plt.show() + # Separate the optimized stream function again onto the different mesh parts for part_ind in range(len(coil_parts)): coil_part = coil_parts[part_ind] diff --git a/pyCoilGen/sub_functions/write_stl_without_sweep.py b/pyCoilGen/sub_functions/write_stl_without_sweep.py new file mode 100644 index 00000000..3e25c04a --- /dev/null +++ b/pyCoilGen/sub_functions/write_stl_without_sweep.py @@ -0,0 +1,257 @@ +# This file ingests the STL file from pyCoilGen that works well only with a sweep +""" +Module: write_stl_without_sweep.py +This module processes STL files generated by pyCoilGen to remove sweep effects from wire paths +and extract centerline coordinates for 3D printing and wire bending applications. +Dependencies: + - trimesh: For loading and manipulating STL mesh files + - numpy: For numerical array operations + - pandas: For data manipulation and Excel file operations + - scipy: For spatial operations and cubic spline interpolation + - scikit-learn: For clustering algorithms (DBSCAN) + - matplotlib: For 3D visualization of wire paths +Required Import for CubicSpline: +Main Workflow: + 1. Load STL file and visualize wire paths with sweep effect + 2. Calculate moving average of x,y coordinates to smooth vertices + 3. Identify deviated coordinates from initial y-plane + 4. Segment curves based on index gaps + 5. Generate smooth spline curves between curve endpoints + 6. Remove duplicate or too-close adjacent vertices + 7. Extract wire path centerline coordinates + 8. Export coordinates to text file for CAD software import +Output: + - Visualization plots of original, adjusted, and smoothed wire paths + - Text file containing wire path coordinates (x, y, z per line) +""" +# and the xyz coordinates of the centers of the wire paths to an excel file. +# Step 1: Load the STL file and visualize the wire paths with sweep +# Step 2: Find the moving average of the x,y coordinates for every 10 vertices +# Step 3: Save the modified mesh to a new STL file +# Step 4: Extract wire path coordinates and save to Excel + +import trimesh +import numpy as np +import pandas as pd +import os +import logging +logging.basicConfig(level=logging.INFO) +log = logging.getLogger(__name__) +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from scipy.spatial import cKDTree +from sklearn.cluster import DBSCAN +from scipy.spatial.distance import cdist +from scipy.interpolate import CubicSpline + +# Step 1: Load the STL file and visualize the wire paths +stl_file_path = "./images/Tenacity_gradient_print_z_wire_1_z.stl" +txt_output_path = stl_file_path.replace('.stl', '.txt') +mesh = trimesh.load_mesh(stl_file_path) +log.info(f"Loaded mesh with {len(mesh.vertices)} vertices and {len(mesh.faces)} faces.") + + +# Now lets visualize the mesh to see the wire paths with sweep from the vertices and faces +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +ax.plot(mesh.vertices[:, 0], mesh.vertices[:, 1], mesh.vertices[:, 2], 'b.') +ax.set_title('With Sweep Wire Paths') +ax.set_xlabel('X axis') +ax.set_ylabel('Y axis') +ax.set_zlabel('Z axis') +plt.show() + + +# Step 2: Find the moving average of the x,y coordinates for every 10 vertices along z to remove the sweep effect +# Calculate the moving average of x, y coordinates for every 10 vertices along z +window_size = 50 +x_coords = mesh.vertices[:, 0] +y_coords = mesh.vertices[:, 1] +z_coords = mesh.vertices[:, 2] +initial_y_plane = np.mean(y_coords[:5]) +vertices_adjusted = mesh.vertices.copy() + +# Create an array to hold the adjusted coordinates +adjusted_coords = vertices_adjusted.copy() +for i in range(len(mesh.vertices)): + start_idx = max(0, i - window_size // 2) + end_idx = min(len(mesh.vertices), i + window_size // 2 + 1) + window = mesh.vertices[start_idx:end_idx, :] + adjusted_coords[i, :] = np.mean(window[:, 0]) + adjusted_coords[i, 1] = np.mean(window[:, 1]) + adjusted_coords[i, 2] = np.mean(window[:, 2]) + + +# The vertices are stacked according to z-levels now but we need them to be back in original like wire paths +vertices_adjusted = adjusted_coords + +# Update mesh with adjusted vertices before checking and exporting +mesh.vertices = vertices_adjusted + + +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +ax.plot(vertices_adjusted[:, 0], vertices_adjusted[:, 1], vertices_adjusted[:, 2], 'b-') +ax.set_title('Adjusted Vertices without Sweep') +ax.set_xlabel('X axis') +ax.set_ylabel('Y axis') +ax.set_zlabel('Z axis') +plt.show() + + +# Identify and plot coordinates in vertices_adjusted that are more than abs(1.75) away from the initial y plane +deviated_coords = vertices_adjusted[np.abs(vertices_adjusted[:, 1] - initial_y_plane) > 0.01] +# Store the corresponding indices of the deviated coordinates +deviated_indices = np.where(np.abs(vertices_adjusted[:, 1] - initial_y_plane) > 0.01)[0] +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +ax.plot(deviated_coords[:, 0], deviated_coords[:, 1], deviated_coords[:, 2], 'ro-') +ax.set_title('Deviated Coordinates from Initial Y Plane') +ax.set_xlabel('X axis') +ax.set_ylabel('Y axis') +ax.set_zlabel('Z axis') +plt.show() + + +# Identify the curves in the deviated coordinates if there is a significant gap between the indices in adjusted vertices, it probably indicates a new curve +# Find gaps in deviated indices to identify separate curves +gaps = np.where(np.diff(deviated_indices) > 1)[0] +curve_starts = np.insert(gaps + 1, 0, 0) # Start of each curve +curve_ends = np.append(gaps, len(deviated_indices) - 1) # End of each curve + +# Log the identified curves +for ind_curve in range(len(curve_starts)): + log.info(f"Identified curve from index {deviated_indices[curve_starts[ind_curve]]} to {deviated_indices[curve_ends[ind_curve]]}") + + +# Change the depth of the three curves to be at different y-plane values to help 3D printing and wire bending +# First smooth the three curves to remove sharp turns in y-axis +normal_shift_length = 2 +curves = [] +sign_use = -1 + +# Plot the smoothed curves +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +for ind_curve in range(len(curve_starts)): + curve_indices = deviated_indices[curve_starts[ind_curve]:curve_ends[ind_curve]+1] + curve_coords = vertices_adjusted[curve_indices, :] + curves.append(curve_coords) + + curve_start = curve_coords[0, :] + curve_end = curve_coords[-1, :] + curve_deepest_y = initial_y_plane + ((normal_shift_length * (ind_curve + 1)) * sign_use) # Shift each curve deeper by 0.5m increments + sign_use *= -1 # Alternate the sign for next curve + + # Generate a smooth spline connecting curve_start and curve_end with a midpoint at curve_deepest_y + t = np.linspace(0, 1, len(curve_coords)) + midpoint = np.array([(curve_start[0] + curve_end[0]) / 2, curve_deepest_y, (curve_start[2] + curve_end[2]) / 2]) + spline_points = np.array([curve_start, midpoint, curve_end]) + + # Create a cubic spline interpolation + cs = CubicSpline([0, 0.5, 1], spline_points, axis=0) + smooth_curve = cs(t) + + # # Ensure the current curve is at least distance 2 away from previous curves + # if ind_curve > 0: + # for prev_curve in curves[:-1]: + # min_dist_to_prev = np.min(cdist(smooth_curve, prev_curve)) + # if min_dist_to_prev < 2: + # # Adjust curve_deepest_y to increase separation + # adjustment = 2 - min_dist_to_prev + 0.5 + # curve_deepest_y += adjustment + # midpoint = np.array([(curve_start[0] + curve_end[0]) / 2, curve_deepest_y, (curve_start[2] + curve_end[2]) / 2]) + # spline_points = np.array([curve_start, midpoint, curve_end]) + # cs = CubicSpline([0, 0.5, 1], spline_points, axis=0) + # smooth_curve = cs(t) + + # Replace the original curve coordinates with the smooth curve + vertices_adjusted[curve_indices, :] = smooth_curve + ax.plot(smooth_curve[:, 0], smooth_curve[:, 1], smooth_curve[:, 2], color=plt.cm.tab10(ind_curve)) + + + +# Sharp turns between curves cause sweep to fail in solidworks - need to smooth those too +# Make sure there are no steep turns > 45 degrees in the x-z plane between adjacent points +angle_store = np.zeros(len(vertices_adjusted)) +for i in range(1, len(vertices_adjusted)-1): + vec1 = vertices_adjusted[i] - vertices_adjusted[i-1] + vec2 = vertices_adjusted[i+1] - vertices_adjusted[i] + cos_angle = np.dot(vec1, vec2) / (np.linalg.norm(vec1) * np.linalg.norm(vec2)) + angle_store[i] = np.rad2deg(np.arccos(np.clip(cos_angle, -1.0, 1.0))) # Angle in radians + +# Plot a 3D figure of the angles at each point to visualize +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +scatter = ax.scatter(vertices_adjusted[:, 0], vertices_adjusted[:, 1], vertices_adjusted[:, 2], c=angle_store, cmap='viridis', s=10) +ax.set_title('Angles at Each Vertex') +ax.set_xlabel('X axis') +ax.set_ylabel('Y axis') +ax.set_zlabel('Z axis') +cbar = plt.colorbar(scatter, ax=ax, label='Angle (degrees)') +plt.show() + + +# Remove the sharp turns by fitting a cubic spline to those segments +sharp_turn_indices = np.where(angle_store > 30)[0] +for idx in sharp_turn_indices: + if idx > 1 and idx < len(vertices_adjusted) - 2: + segment_indices = np.arange(idx - 2, idx + 3) + segment_coords = vertices_adjusted[segment_indices, :] + t_segment = np.linspace(0, 1, len(segment_coords)) + cs_segment = CubicSpline(t_segment, segment_coords, axis=0) + t_fine = np.linspace(0, 1, 50) + smooth_segment = cs_segment(t_fine) + # Replace the original segment with the smooth segment (interpolated back to original indices) + for j, seg_idx in enumerate(segment_indices): + vertices_adjusted[seg_idx, :] = smooth_segment[int(j * (len(t_fine) - 1) / (len(segment_indices) - 1)), :] + +# Make sure the vertices are all 2 decimal places before writing to the text file +vertices_adjusted = np.round(vertices_adjusted, 2) + +# Step 4: Remove duplicate or too-close adjacent vertices +min_distance = 2 # Minimum distance threshold +filtered_vertices = [vertices_adjusted[0]] + +for i in range(1, len(vertices_adjusted)): + distance = np.linalg.norm(vertices_adjusted[i] - filtered_vertices[-1]) + if distance > min_distance: + filtered_vertices.append(vertices_adjusted[i]) + +vertices_adjusted = np.array(filtered_vertices) +log.info(f"Removed {len(filtered_vertices) - len(vertices_adjusted)} duplicate or too-close vertices") + + + +print(f"Number of vertices after interpolation: {len(vertices_adjusted)}") +# Interpolate to get new vertices +# Visualize the smoothed wire path +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +ax.plot(vertices_adjusted[:, 0], vertices_adjusted[:, 1], vertices_adjusted[:, 2], 'r-') +ax.set_title('Final Smoothed Wire Path without Sweep') +ax.set_xlabel('X axis') +ax.set_ylabel('Y axis') +ax.set_zlabel('Z axis') +plt.show() + + + + +# Extract wire path centerline coordinates +wire_path = np.array([[np.round(vertices_adjusted[i, 0], 2), np.round(vertices_adjusted[i, 1], 2), np.round(vertices_adjusted[i, 2], 2)] + for i in range(len(vertices_adjusted)-5)]) # will skip the last vertex to enable operations in solidworks + + +# Also write a simple text file for direct import to solidworks with no heading but only coordinates - one coordianate per line and x, y and z separated by space + +with open(txt_output_path, 'w') as f: + for point in wire_path: + f.write(f"{point[0]} {point[1]} {point[2]}\n") +log.info(f"Saved wire path coordinates to {txt_output_path}") + + + + + + diff --git a/requirements.txt b/requirements.txt index c2f68230..df9432d7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,18 +1,93 @@ -numpy==1.* - # Note: Also need BLAS and gfortran to install scipy # sudo apt-get install libopenblas-dev # sudo apt install gfortran scipy==1.* # For Mesh support -trimesh==3.* -# Missing Trimesh dependency + networkx==2.* # For target field calculations sympy==1.* # For visualisation of matrices -pillow<=9.5 -matplotlib==3.* \ No newline at end of file +matplotlib==3.* +certifi==2026.2.25 +charset-normalizer==3.4.6 +docopt==0.6.2 +idna==3.11 +pipreqs==0.4.13 +requests==2.32.5 +urllib3==2.6.3 +yarg==0.1.10 +attrs==26.1.0 +certifi==2026.2.25 +charset-normalizer==3.4.6 +contourpy==1.3.3 +cycler==0.12.1 +cyclopts==4.10.1 +docopt==0.6.2 +docstring_parser==0.17.0 +docutils==0.22.4 +fonttools==4.62.1 +idna==3.11 +kiwisolver==1.5.0 +markdown-it-py==4.0.0 +matplotlib==3.10.8 +mdurl==0.1.2 +numpy==2.4.3 +packaging==26.0 +pillow==12.1.1 +pipreqs==0.4.13 +platformdirs==4.9.4 +pooch==1.9.0 +Pygments==2.19.2 +pyparsing==3.3.2 +python-dateutil==2.9.0.post0 +pyvista==0.47.1 +requests==2.32.5 +rich==14.3.3 +rich-rst==1.3.2 +scooby==0.11.0 +six==1.17.0 +typing_extensions==4.15.0 +urllib3==2.6.3 +vtk==9.6.0 +yarg==0.1.10 +attrs==26.1.0 +certifi==2026.2.25 +charset-normalizer==3.4.6 +contourpy==1.3.3 +cycler==0.12.1 +cyclopts==4.10.1 +docopt==0.6.2 +docstring_parser==0.17.0 +docutils==0.22.4 +fonttools==4.62.1 +idna==3.11 +kiwisolver==1.5.0 +markdown-it-py==4.0.0 +matplotlib==3.10.8 +mdurl==0.1.2 +numpy==2.4.3 +packaging==26.0 +pillow==12.1.1 +pipreqs==0.4.13 +platformdirs==4.9.4 +pooch==1.9.0 +Pygments==2.19.2 +pyparsing==3.3.2 +python-dateutil==2.9.0.post0 +pyvista==0.47.1 +requests==2.32.5 +rich==14.3.3 +rich-rst==1.3.2 +scooby==0.11.0 +six==1.17.0 +typing_extensions==4.15.0 +urllib3==2.6.3 +vtk==9.6.0 +yarg==0.1.10 +magpylib==5.2.2 +trimesh==4.11.4 +manifold3d==3.4.0 diff --git a/tests/test_stl_generation.py b/tests/test_stl_generation.py new file mode 100644 index 00000000..1dc8463a --- /dev/null +++ b/tests/test_stl_generation.py @@ -0,0 +1,67 @@ +import numpy as np +import trimesh +from scipy.interpolate import splprep, splev + +# ------------------------ +# LOAD + SMOOTH PATH +# ------------------------ +data = np.loadtxt('images/biplanar_xgradient_Tenacity_x_plate0_loop0.txt') + +x, y = data[:, 0], data[:, 1] + +def smooth_path(x, y, smoothing=5.0, n_points=300): + tck, _ = splprep([x, y], s=smoothing, per=True) + u = np.linspace(0, 1, n_points) + x_new, y_new = splev(u, tck) + return np.array(x_new), np.array(y_new) + +x, y = smooth_path(x, y) + +points = np.column_stack([x, y, np.zeros_like(x)]) + +# ------------------------ +# CREATE PLATE +# ------------------------ +plate = trimesh.creation.cylinder( + radius=76.2, + height=3, + sections=128 +) + +plate.apply_translation([0, 0, 1.5]) + +# ------------------------ +# CREATE TUBE (ROBUST) +# ------------------------ +radius = 1.6 +segments = [] + +for i in range(len(points) - 1): + p1 = points[i] + p2 = points[i + 1] + + segment = trimesh.creation.cylinder( + radius=radius, + segment=[p1, p2] + ) + + segments.append(segment) + +tube = trimesh.util.concatenate(segments) + +# Embed into plate +tube.apply_translation([0, 0, 2]) + +# ------------------------ +# BOOLEAN DIFFERENCE +# ------------------------ +result = plate.difference(tube) + +# ------------------------ +# EXPORT +# ------------------------ +result.export('plate_trimesh.stl') + +print("✅ STL generated with Trimesh") + +