Skip to content

Instantly share code, notes, and snippets.

@brossi
Last active April 9, 2025 16:13
Show Gist options
  • Save brossi/e1ff63cc2e4537c4d38b58609036c959 to your computer and use it in GitHub Desktop.
Save brossi/e1ff63cc2e4537c4d38b58609036c959 to your computer and use it in GitHub Desktop.
Material Spec Library
"""Utilities for the Radiolucent Material and Marker Optimization Tool."""
import numpy as np
import math
from typing import List, Dict, Tuple, Union, Optional
def calculate_marker_volume(diameter: float) -> float:
"""Calculate the volume of a spherical marker in mm³
Args:
diameter: The marker diameter in mm
Returns:
The volume in mm³
"""
return (4.0 / 3.0) * np.pi * (diameter / 2.0)**3
def calculate_total_marker_volume(markers: List) -> Tuple[float, float]:
"""Calculate the total volume of all markers in mm³ and cm³
Args:
markers: List of Marker objects
Returns:
A tuple of (volume_mm3, volume_cm3)
"""
if not markers:
return 0.0, 0.0
volume_mm3 = sum([calculate_marker_volume(m.diameter) for m in markers])
volume_cm3 = volume_mm3 / 1000.0
return volume_mm3, volume_cm3
def calculate_panel_weight(panel_length: float, panel_width: float, panel_thickness: float,
substrate_density: float, markers: List = None,
marker_material_density: Optional[float] = None) -> float:
"""Calculate the weight of a panel in grams, including markers if provided
Args:
panel_length: Panel length in mm
panel_width: Panel width in mm
panel_thickness: Panel thickness in mm
substrate_density: Substrate material density in g/cm³
markers: Optional list of Marker objects
marker_material_density: Density of marker material in g/cm³, required if markers are provided
Returns:
Weight in grams
"""
# Convert dimensions from mm to cm
volume_cm3 = (panel_length/10.0) * (panel_width/10.0) * (panel_thickness/10.0)
# Base substrate weight
weight_g = volume_cm3 * substrate_density
# Adjust for markers if present
if markers and marker_material_density is not None:
_, total_marker_volume_cm3 = calculate_total_marker_volume(markers)
# Subtract weight of displaced substrate, add weight of markers
weight_g -= total_marker_volume_cm3 * substrate_density
weight_g += total_marker_volume_cm3 * marker_material_density
return weight_g
def find_max_stress_concentration(stress_results: Dict) -> float:
"""Find the maximum finite stress concentration factor from results
Args:
stress_results: Dictionary of stress concentration calculation results
Returns:
The maximum finite stress concentration factor, or 1.0 if none found
"""
if not stress_results:
return 1.0
finite_kts = [res['stress_concentration_factor']
for res in stress_results.values()
if res['stress_concentration_factor'] != float('inf')]
if finite_kts:
return max(finite_kts)
elif stress_results: # All results have infinite Kt
return float('inf')
else:
return 1.0
def check_material_keys(materials_dict: Dict,
substrate_key: str = None, marker_key: str = None) -> Tuple[bool, str]:
"""Check if material keys are valid and appropriate for their roles
Args:
materials_dict: Dictionary of materials
substrate_key: Key for the substrate material (should be radiolucent)
marker_key: Key for the marker material (should be non-radiolucent)
Returns:
Tuple of (is_valid, error_message)
"""
if substrate_key and substrate_key not in materials_dict:
return False, f"Substrate material key '{substrate_key}' not found in library."
if marker_key and marker_key not in materials_dict:
return False, f"Marker material key '{marker_key}' not found in library."
if substrate_key and not materials_dict[substrate_key].is_radiolucent:
return False, f"Selected substrate material '{substrate_key}' is not radiolucent."
if marker_key and materials_dict[marker_key].is_radiolucent:
return False, f"Selected marker material '{marker_key}' is radiolucent."
return True, ""
def calculate_linear_attenuation(material) -> float:
"""Calculate the linear attenuation coefficient of a material
Args:
material: A Material object with x_ray_attenuation and density properties
Returns:
Linear attenuation coefficient in cm^-1
"""
return material.x_ray_attenuation * material.density
def filter_materials_by_type(materials_dict: Dict, is_radiolucent: bool) -> Dict:
"""Filter materials by their radiolucency
Args:
materials_dict: Dictionary of Material objects
is_radiolucent: True to get substrate materials, False for marker materials
Returns:
Filtered dictionary of materials
"""
return {k: v for k, v in materials_dict.items() if v.is_radiolucent == is_radiolucent}
import pandas as pd
import numpy as np
from typing import Dict, List, Tuple, Optional
from utils.panel import SubstratePanel
from utils.materials import Material, Marker
from utils.calculations import (
calculate_panel_weight,
find_max_stress_concentration,
check_material_keys,
filter_materials_by_type
)
def compare_substrate_materials(materials_dict: Dict[str, Material],
marker_material_key: str,
marker_diameter: float,
panel_dimensions: Tuple[float, float],
load: float,
num_markers: int,
safety_factor: float = 2.0) -> pd.DataFrame:
"""Compares different substrate materials for a fixed marker configuration.
Calculates the minimum required thickness for each substrate based on load
and the simplified stress model, then evaluates radiolucency and weight.
Args:
materials_dict: Dictionary of all available Material objects.
marker_material_key: Key (str) of the marker material to use from materials_dict.
marker_diameter: Diameter of markers in mm.
panel_dimensions: Tuple (length, width) of the panel in mm.
load: Total load on the panel in Newtons (N).
num_markers: Number of markers to place in a grid pattern.
safety_factor: Safety factor to use for thickness calculation.
Returns:
pandas DataFrame comparing substrate materials.
"""
results = []
# Check for valid marker material
is_valid, error_msg = check_material_keys(materials_dict, marker_key=marker_material_key)
if not is_valid:
print(f"Error: {error_msg}")
return pd.DataFrame()
marker_material = materials_dict[marker_material_key]
# Filter for substrate materials only
substrate_materials = filter_materials_by_type(materials_dict, is_radiolucent=True)
if not substrate_materials:
print("Error: No radiolucent substrate materials found in library.")
return pd.DataFrame()
print(f"Comparing {len(substrate_materials)} substrate materials...")
for mat_id, material in substrate_materials.items():
print(f" Analyzing substrate: {material.name}")
# Create substrate with a nominal starting thickness (will be recalculated)
panel = SubstratePanel(
material=material,
length=panel_dimensions[0],
width=panel_dimensions[1],
thickness=marker_diameter * 1.1 # Initial guess > marker diameter
)
# Generate marker grid pattern
markers = panel.generate_grid_pattern(
marker_material=marker_material,
marker_diameter=marker_diameter,
num_markers=num_markers
)
panel.markers = markers
# Calculate minimum thickness required for this substrate
min_thickness = panel.calculate_minimum_thickness(load, safety_factor=safety_factor)
# Update panel with minimum required thickness and recalculate properties
if min_thickness == float('inf'):
print(f" Skipping {material.name}: Infinite thickness required (likely marker overlap).")
# Add result with NaN or Inf values
result = {
'Substrate': material.name,
'Min Thickness (mm)': np.inf,
'Radiolucency Score': np.nan,
'Weight (g)': np.inf,
'Flexural Strength (MPa)': material.flexural_strength,
'Young\'s Modulus (GPa)': material.youngs_modulus
}
else:
panel.thickness = min_thickness
radiolucency = panel.calculate_radiolucency()
# Calculate weight
weight_g = calculate_panel_weight(
panel_length=panel.length,
panel_width=panel.width,
panel_thickness=panel.thickness,
substrate_density=material.density,
markers=panel.markers,
marker_material_density=marker_material.density
)
result = {
'Substrate': material.name,
'Min Thickness (mm)': round(min_thickness, 3),
'Radiolucency Score': round(radiolucency, 2),
'Weight (g)': round(weight_g, 2),
'Flexural Strength (MPa)': material.flexural_strength,
'Young\'s Modulus (GPa)': material.youngs_modulus
}
print(f" Min Thickness: {min_thickness:.2f} mm, Radiolucency: {radiolucency:.1f}, Weight: {weight_g:.1f} g")
results.append(result)
if not results:
return pd.DataFrame() # Return empty if no substrates could be analyzed
df = pd.DataFrame(results)
df = df.sort_values(by='Radiolucency Score', ascending=False).reset_index(drop=True)
return df
def compare_marker_materials(substrate_material_key: str,
materials_dict: Dict[str, Material],
marker_diameter: float,
panel_dimensions: Tuple[float, float],
thickness: float,
num_markers: int) -> pd.DataFrame:
"""Compares different marker materials within a fixed substrate and geometry.
Evaluates radiolucency, relative visibility, and the simplified max stress
concentration factor for each marker material.
Args:
substrate_material_key: Key (str) of the substrate material from materials_dict.
materials_dict: Dictionary of all available Material objects.
marker_diameter: Diameter of markers in mm.
panel_dimensions: Tuple (length, width) of the panel in mm.
thickness: The *fixed* thickness of the panel in mm for this comparison.
num_markers: Number of markers to place in a grid pattern.
Returns:
pandas DataFrame comparing marker materials.
"""
results = []
# Check for valid substrate material
is_valid, error_msg = check_material_keys(materials_dict, substrate_key=substrate_material_key)
if not is_valid:
print(f"Error: {error_msg}")
return pd.DataFrame()
substrate_material = materials_dict[substrate_material_key]
# Filter for non-radiolucent (marker) materials only
marker_materials = filter_materials_by_type(materials_dict, is_radiolucent=False)
if not marker_materials:
print("Error: No non-radiolucent marker materials found in library.")
return pd.DataFrame()
# Create base panel
panel = SubstratePanel(
material=substrate_material,
length=panel_dimensions[0],
width=panel_dimensions[1],
thickness=thickness
)
# Generate one marker pattern (positions) using a temporary marker material
# The positions will be the same for all marker materials being compared.
temp_marker_mat = list(marker_materials.values())[0]
base_markers_positions = [m.position for m in panel.generate_grid_pattern(
marker_material=temp_marker_mat,
marker_diameter=marker_diameter,
num_markers=num_markers
)]
print(f"Comparing {len(marker_materials)} marker materials...")
for mat_id, material in marker_materials.items():
# Create new list of markers with the current material but same positions
current_markers = [Marker(material=material, diameter=marker_diameter, position=pos)
for pos in base_markers_positions]
# Update panel with these markers
panel.markers = current_markers
# Calculate properties
radiolucency = panel.calculate_radiolucency()
stress_results = panel.calculate_stress_concentration()
max_kt = find_max_stress_concentration(stress_results)
# Relative Visibility: Higher density and attenuation coefficient mean more visible
# Linear attenuation coeff (mu) = mass_attenuation * density
relative_visibility = material.x_ray_attenuation * material.density
result = {
'Marker Material': material.name,
'Lin. Attenuation (cm⁻¹)': round(relative_visibility, 3),
'Density (g/cm³)': material.density,
'Radiolucency Score': round(radiolucency, 2),
'Max Stress Conc. (Kt)': round(max_kt, 2) if max_kt != float('inf') else np.inf,
'Relative Visibility Score': round(relative_visibility * 10, 2) # Scale for better comparison range
}
print(f" Analyzed marker: {material.name}, Radiolucency: {radiolucency:.1f}, Visibility Score: {result['Relative Visibility Score']:.1f}")
results.append(result)
if not results:
return pd.DataFrame()
df = pd.DataFrame(results)
# Sort by a balance of visibility and radiolucency? Maybe let user decide.
df = df.sort_values(by='Radiolucency Score', ascending=False).reset_index(drop=True)
return df
def design_and_analyze_panel(substrate_material_key: str,
marker_material_key: str,
marker_diameter: float,
panel_length: float,
panel_width: float,
num_markers: int,
load: float,
safety_factor: float = 2.0,
min_edge_distance: float = 5.0,
materials_dict: Dict = None) -> Optional[SubstratePanel]:
"""Creates a panel with specified parameters, analyzes it, and visualizes.
Calculates minimum thickness, radiolucency, weight, and shows the layout
with simplified stress concentration visualization.
Args:
substrate_material_key: Key of the substrate material.
marker_material_key: Key of the marker material.
marker_diameter: Diameter of markers in mm.
panel_length: Panel length in mm.
panel_width: Panel width in mm.
num_markers: Number of markers for the grid pattern.
load: Total load in Newtons (N).
safety_factor: Safety factor for thickness calculation.
min_edge_distance: Minimum distance from marker center to edge (mm).
materials_dict: The dictionary containing material objects.
Returns:
The analyzed SubstratePanel object, or None if inputs are invalid.
"""
if materials_dict is None:
print("Error: No materials dictionary provided.")
return None
# Check materials keys
is_valid, error_msg = check_material_keys(
materials_dict,
substrate_key=substrate_material_key,
marker_key=marker_material_key
)
if not is_valid:
print(f"Error: {error_msg}")
return None
substrate_mat = materials_dict[substrate_material_key]
marker_mat = materials_dict[marker_material_key]
# Create the panel with a nominal thickness first
panel = SubstratePanel(
material=substrate_mat,
length=panel_length,
width=panel_width,
thickness=marker_diameter * 1.1 # Initial guess > marker diameter
)
# Generate grid pattern markers
markers = panel.generate_grid_pattern(
marker_material=marker_mat,
marker_diameter=marker_diameter,
num_markers=num_markers,
min_edge_distance=min_edge_distance
)
panel.markers = markers
# Calculate minimum thickness
min_thickness = panel.calculate_minimum_thickness(load, safety_factor)
if min_thickness == float('inf'):
print("Error: Cannot determine finite minimum thickness. Check marker placement/overlap.")
# Visualize the problematic layout anyway
panel.thickness = marker_diameter * 1.1 # Use guess thickness for viz
fig, ax = panel.visualize(show_stress=True)
plt.suptitle("Panel Layout (Infinite Thickness Required)", y=1.02)
plt.show()
return panel # Return the panel object even if thickness calc failed
panel.thickness = min_thickness
# Calculate other properties with final thickness
radiolucency = panel.calculate_radiolucency()
# Weight calculation
weight_g = calculate_panel_weight(
panel_length=panel.length,
panel_width=panel.width,
panel_thickness=panel.thickness,
substrate_density=panel.material.density,
markers=panel.markers,
marker_material_density=marker_mat.density
)
# Display results
print("--- Panel Design Analysis ---")
print(f"Substrate: {panel.material.name} ({substrate_material_key})")
print(f"Marker Material: {marker_mat.name} ({marker_material_key})")
print(f"Panel Size: {panel_length} L x {panel_width} W mm")
print(f"Markers: {num_markers} x {marker_diameter}mm diameter (Grid Pattern)")
print(f"Load: {load} N, Safety Factor: {safety_factor}")
print("-----------------------------")
print(f"Est. Minimum Thickness: {min_thickness:.3f} mm")
print(f"Radiolucency Score: {radiolucency:.2f}")
print(f"Estimated Weight: {weight_g:.2f} g")
print("-----------------------------")
print("Note: Stress visualization uses a simplified model.")
# Visualize
import matplotlib.pyplot as plt
fig, ax = panel.visualize(show_stress=True)
plt.suptitle("Panel Layout and Simplified Stress Concentration", y=1.02)
plt.show()
return panel
import requests
import json
from typing import Dict, List, Optional
from dataclasses import dataclass
from typing import Tuple
@dataclass
class Material:
"""Represents substrate and marker materials.
Attributes:
name: Name of the material.
density: Material density in g/cm³.
x_ray_attenuation: Mass attenuation coefficient in cm²/g (typically at a specific energy like 100 keV).
youngs_modulus: Young's Modulus (stiffness) in GPa.
flexural_strength: Flexural strength (bending strength) in MPa.
is_radiolucent: Boolean indicating if the material is considered radiolucent (True for substrates).
"""
name: str
density: float # g/cm³
x_ray_attenuation: float # cm²/g at 100 keV (assumed)
youngs_modulus: float # GPa
flexural_strength: float # MPa
is_radiolucent: bool # True for substrate materials, False for markers
@dataclass
class Marker:
"""Represents a single spherical marker.
Attributes:
material: The Material object defining the marker's properties.
diameter: The diameter of the spherical marker in mm.
position: A tuple representing the (x, y) coordinates of the marker's center in mm.
"""
material: Material
diameter: float # mm
position: Tuple[float, float] # (x, y) center position in mm
def load_materials_from_url(url: str) -> Dict[str, Material]:
"""Loads material data from a JSON file at the specified URL.
Args:
url: The URL to the JSON file containing material data.
Returns:
A dictionary mapping material keys (str) to Material objects.
Returns an empty dictionary if loading fails.
"""
materials_dict = {}
try:
response = requests.get(url)
response.raise_for_status() # Raise an exception for bad status codes (4xx or 5xx)
data = json.loads(response.text)
for key, props in data.items():
try:
materials_dict[key] = Material(
name=props['name'],
density=props['density'],
x_ray_attenuation=props['x_ray_attenuation'],
youngs_modulus=props['youngs_modulus'],
flexural_strength=props['flexural_strength'],
is_radiolucent=props['is_radiolucent']
)
except KeyError as e:
print(f"Warning: Missing key {e} for material '{key}'. Skipping.")
except Exception as e:
print(f"Warning: Error processing material '{key}': {e}. Skipping.")
print(f"Successfully loaded {len(materials_dict)} materials from URL.")
except requests.exceptions.RequestException as e:
print(f"Error fetching materials from URL: {e}")
print("Please check the URL and your internet connection.")
return {}
except json.JSONDecodeError as e:
print(f"Error decoding JSON from URL: {e}")
return {}
except Exception as e:
print(f"An unexpected error occurred during material loading: {e}")
return {}
return materials_dict
def display_available_materials(materials_dict: Dict[str, Material]) -> None:
"""Display the available materials separated by type
Args:
materials_dict: Dictionary of Material objects
"""
if not materials_dict:
print("\nMaterial library could not be loaded or is empty.")
return
print("\nAvailable Substrate Materials:")
for mat_id, material in materials_dict.items():
if material.is_radiolucent:
print(f"- {material.name} (Key: {mat_id})")
print("\nAvailable Marker Materials:")
for mat_id, material in materials_dict.items():
if not material.is_radiolucent:
print(f"- {material.name} (Key: {mat_id})")
def get_default_materials() -> Dict[str, Material]:
"""Create a basic set of default materials if loading from URL fails
Returns:
Dictionary of default Material objects
"""
return {
# Substrate materials (radiolucent)
"carbon_fiber_composite": Material(
name="Carbon Fiber Composite",
density=1.6,
x_ray_attenuation=0.18,
youngs_modulus=70.0,
flexural_strength=600.0,
is_radiolucent=True
),
"acrylic": Material(
name="Acrylic (PMMA)",
density=1.18,
x_ray_attenuation=0.20,
youngs_modulus=3.0,
flexural_strength=100.0,
is_radiolucent=True
),
# Marker materials (radiopaque)
"ss304": Material(
name="Stainless Steel 304",
density=8.0,
x_ray_attenuation=2.0,
youngs_modulus=193.0,
flexural_strength=860.0,
is_radiolucent=False
),
"tungsten": Material(
name="Tungsten",
density=19.3,
x_ray_attenuation=3.5,
youngs_modulus=400.0,
flexural_strength=1510.0,
is_radiolucent=False
)
}
{
"pei_30gf": {
"name": "PEI with 30% Glass Fiber",
"density": 1.51,
"x_ray_attenuation": 0.22,
"youngs_modulus": 9.3,
"flexural_strength": 170,
"is_radiolucent": true
},
"pei_cf": {
"name": "PEI with Carbon Fiber",
"density": 1.41,
"x_ray_attenuation": 0.17,
"youngs_modulus": 22.0,
"flexural_strength": 300,
"is_radiolucent": true
},
"carbon_fiber_composite": {
"name": "Carbon Fiber Composite",
"density": 1.6,
"x_ray_attenuation": 0.18,
"youngs_modulus": 70.0,
"flexural_strength": 600.0,
"is_radiolucent": true
},
"peek_cf": {
"name": "PEEK with Carbon Fiber",
"density": 1.32,
"x_ray_attenuation": 0.16,
"youngs_modulus": 18.0,
"flexural_strength": 240,
"is_radiolucent": true
},
"acrylic": {
"name": "Acrylic (PMMA)",
"density": 1.18,
"x_ray_attenuation": 0.20,
"youngs_modulus": 3.0,
"flexural_strength": 100.0,
"is_radiolucent": true
},
"ss304": {
"name": "Stainless Steel 304",
"density": 8.0,
"x_ray_attenuation": 2.0,
"youngs_modulus": 193.0,
"flexural_strength": 860.0,
"is_radiolucent": false
},
"copper": {
"name": "Copper",
"density": 8.96,
"x_ray_attenuation": 1.6,
"youngs_modulus": 110.0,
"flexural_strength": 330.0,
"is_radiolucent": false
},
"zirconia": {
"name": "Zirconia Ceramic",
"density": 6.05,
"x_ray_attenuation": 0.9,
"youngs_modulus": 200.0,
"flexural_strength": 900.0,
"is_radiolucent": false
},
"tungsten": {
"name": "Tungsten",
"density": 19.3,
"x_ray_attenuation": 3.5,
"youngs_modulus": 400.0,
"flexural_strength": 1510.0,
"is_radiolucent": false
}
}
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import math
from typing import List, Dict, Tuple, Optional
from utils.materials import Material, Marker
from utils.calculations import calculate_marker_volume, calculate_total_marker_volume
class SubstratePanel:
"""Represents a substrate panel with embedded spherical markers.
Handles calculations for radiolucency, simplified stress concentration,
minimum thickness estimation, marker pattern generation, and visualization.
Attributes:
material: The Material object for the substrate.
length: The length of the panel (x-dimension) in mm.
width: The width of the panel (y-dimension) in mm.
thickness: The thickness of the panel (z-dimension) in mm.
markers: A list of Marker objects embedded in the panel.
"""
def __init__(self,
material: Material,
length: float, # mm
width: float, # mm
thickness: float, # mm
markers: List[Marker] = None):
"""Initializes the SubstratePanel.
Args:
material: The substrate material.
length: Panel length in mm.
width: Panel width in mm.
thickness: Panel thickness in mm.
markers: Optional list of markers to initialize with.
"""
if not material.is_radiolucent:
print(f"Warning: Initializing SubstratePanel with non-radiolucent material '{material.name}'.")
self.material = material
self.length = length
self.width = width
self.thickness = thickness
self.markers = markers if markers is not None else []
def add_marker(self, marker: Marker):
"""Adds a single marker to the panel's list of markers."""
if marker.diameter > self.thickness:
print(f"Warning: Marker diameter ({marker.diameter}mm) exceeds panel thickness ({self.thickness}mm).")
self.markers.append(marker)
def calculate_radiolucency(self) -> float:
"""Calculates the overall radiolucency score of the panel.
This score represents the estimated X-ray transmission percentage (0-100)
based on the Beer-Lambert law, considering both substrate and spherical
marker attenuation (assumed at 100 keV).
Higher score means better radiolucency (less attenuation).
Returns:
The calculated radiolucency score (0-100).
"""
# Linear attenuation coefficient (mu) = mass_attenuation * density
# Attenuation = mu * path_length (where path_length is thickness in cm)
substrate_mu = self.material.x_ray_attenuation * self.material.density # cm^-1
substrate_attenuation_factor = substrate_mu * (self.thickness / 10.0) # Unitless exponent
total_panel_volume_mm3 = self.length * self.width * self.thickness
total_marker_volume_mm3, _ = calculate_total_marker_volume(self.markers)
weighted_marker_attenuation_factor = 0
if not self.markers:
# No markers, only substrate attenuation
total_attenuation_factor = substrate_attenuation_factor
else:
for marker in self.markers:
# Volume of a sphere
volume_mm3 = calculate_marker_volume(marker.diameter)
# Marker attenuation factor calculation
marker_mu = marker.material.x_ray_attenuation * marker.material.density # cm^-1
# Assume path length through marker is its diameter for simplicity in avg
# This is an approximation for averaging purposes.
marker_attenuation_factor = marker_mu * (marker.diameter / 10.0) # Unitless exponent
# Weight marker attenuation by its volume fraction contribution
weighted_marker_attenuation_factor += marker_attenuation_factor * (volume_mm3 / total_panel_volume_mm3)
# Calculate substrate volume fraction
substrate_volume_fraction = 1.0 - (total_marker_volume_mm3 / total_panel_volume_mm3)
# Calculate total average attenuation factor (weighted average)
total_attenuation_factor = (substrate_attenuation_factor * substrate_volume_fraction) + weighted_marker_attenuation_factor
# Radiolucency = Transmission percentage = 100 * exp(-total_attenuation_factor)
# Handle potential overflow if attenuation is extremely high
try:
radiolucency_score = 100.0 * math.exp(-total_attenuation_factor)
except OverflowError:
radiolucency_score = 0.0 # Essentially zero transmission
return max(0.0, min(100.0, radiolucency_score)) # Clamp between 0 and 100
def calculate_stress_concentration(self) -> Dict:
"""Calculates simplified stress concentration factors (Kt) near each marker.
*** WARNING: Uses a highly simplified model! ***
The formula Kt = 1 + 2 * (marker.diameter / s) is a basic approximation
and does NOT account for material property mismatch (marker vs. substrate),
3D effects, finite panel dimensions, complex loading, or interactions
between multiple markers accurately. Results should be used for RELATIVE
comparison only. Consult FEA for accurate stress analysis.
Here, 's' is the minimum distance from the marker's surface to the nearest
edge or the surface of the nearest neighboring marker in the 2D plane.
Returns:
A dictionary mapping marker index (int) to details:
{'marker': Marker, 'stress_concentration_factor': float, 'nearest_distance': float}
"""
results = {}
if not self.markers:
return results
positions = np.array([m.position for m in self.markers])
diameters = np.array([m.diameter for m in self.markers])
radii = diameters / 2.0
for i, marker in enumerate(self.markers):
pos_i = positions[i]
radius_i = radii[i]
# Calculate minimum distance to other markers' surfaces
min_marker_distance = float('inf')
if len(self.markers) > 1:
# Calculate distances to all other marker centers
diffs = positions - pos_i
center_distances = np.sqrt(np.sum(diffs**2, axis=1))
# Calculate surface-to-surface distances (ignore self, i.e., distance=0)
surface_distances = center_distances - (radii + radius_i)
# Find minimum distance to *another* marker's surface
valid_distances = surface_distances[surface_distances > 1e-9] # Avoid self-comparison or touching markers giving zero
if len(valid_distances) > 0:
min_marker_distance = np.min(valid_distances)
else: # Handle cases where markers might be overlapping exactly or only one marker exists
# Check if any center distance (other than self) is <= sum of radii
if np.any( (center_distances > 1e-9) & (center_distances <= (radii + radius_i)) ):
min_marker_distance = 0 # Touching or overlapping
# else: min_marker_distance remains inf if only one marker or widely separated
# Calculate minimum distance to edges (center to edge minus radius)
dist_to_left = pos_i[0] - radius_i
dist_to_right = (self.length - pos_i[0]) - radius_i
dist_to_bottom = pos_i[1] - radius_i
dist_to_top = (self.width - pos_i[1]) - radius_i
min_edge_distance = min(dist_to_left, dist_to_right, dist_to_bottom, dist_to_top)
# Overall minimum distance 's'
s = min(min_marker_distance, min_edge_distance)
# Calculate simplified Kt
if s <= 1e-9: # Use tolerance for floating point comparisons
kt = float('inf') # Markers are overlapping or touching edge/each other
s = 0 # Report nearest distance as 0
else:
# Simplified formula - USE WITH CAUTION
kt = 1.0 + 2.0 * (marker.diameter / s)
results[i] = {
'marker': marker,
'stress_concentration_factor': kt,
'nearest_distance': s
}
return results
def calculate_minimum_thickness(self, load: float, safety_factor: float = 2.0) -> float:
"""Estimates the minimum required panel thickness based on load and marker placement.
Assumes the panel is a simply supported plate under uniform pressure.
Uses the simplified stress concentration factor (Kt) calculated previously.
Ensures thickness is sufficient to fully embed the largest marker.
Args:
load: Total load applied uniformly over the panel surface in Newtons (N).
safety_factor: The desired factor of safety against flexural failure.
Returns:
The estimated minimum thickness in mm. Returns float('inf') if markers
are overlapping or touching edges (infinite Kt).
"""
# Get maximum *finite* stress concentration factor
stress_results = self.calculate_stress_concentration()
max_kt = 1.0
finite_kts = [res['stress_concentration_factor']
for res in stress_results.values()
if res['stress_concentration_factor'] != float('inf')]
if finite_kts:
max_kt = max(finite_kts)
elif stress_results: # If there are markers but all have infinite Kt
print("Warning: Markers are overlapping or touching edges. Cannot calculate finite minimum thickness based on stress.")
max_kt = float('inf')
# Proceed to calculate thickness based on marker size only
# --- Thickness based on mechanical strength ---
min_thickness_mech = 0
if max_kt == float('inf'):
min_thickness_mech = float('inf')
else:
# Calculate bending stress for a simply supported rectangular plate under uniform load
# Formula: σ_max ≈ β * q * b² / t² (where β depends on a/b ratio, often ~0.5-0.75 for center)
# Using a simplified form: σ_max ≈ 0.75 * q * b² / t² for conservative estimate
# where q = pressure (N/mm² = MPa), b = shorter dimension (mm), t = thickness (mm)
# Convert load (N) to pressure (MPa = N/mm²)
panel_area_mm2 = self.length * self.width
if panel_area_mm2 <= 0:
print("Warning: Panel area is zero or negative. Cannot calculate thickness.")
return float('inf')
pressure_mpa = load / panel_area_mm2
shorter_dim_mm = min(self.length, self.width)
# Allowable stress = Flexural Strength / (Kt * Safety Factor)
# Note: Flexural strength is usually in MPa, which is N/mm²
allowable_stress_mpa = self.material.flexural_strength / (max_kt * safety_factor)
if allowable_stress_mpa <= 0:
print("Warning: Allowable stress is zero or negative. Cannot calculate thickness.")
# This could happen if flexural strength is zero or Kt is infinite (already handled)
min_thickness_mech = float('inf')
else:
# Solve for minimum thickness t: t² = 0.75 * pressure * b² / allowable_stress
# Coefficient 0.75 is approximate, depends on exact boundary conditions and aspect ratio.
coefficient = 0.75
min_thickness_mech_squared = (coefficient * pressure_mpa * shorter_dim_mm**2) / allowable_stress_mpa
if min_thickness_mech_squared < 0:
print("Warning: Intermediate calculation for thickness squared is negative.")
min_thickness_mech = float('inf')
else:
min_thickness_mech = np.sqrt(min_thickness_mech_squared)
# --- Thickness based on accommodating markers ---
min_thickness_geom = 0
if self.markers:
max_marker_diameter = max(marker.diameter for marker in self.markers)
# Ensure thickness is at least large enough to contain the marker sphere.
# Adding a buffer (e.g., 1.5x) ensures some substrate material above/below.
# The exact requirement depends on manufacturing/embedding process.
min_thickness_geom = max_marker_diameter * 1.1 # Ensure at least diameter + 10% buffer
# Final minimum thickness is the max of the two requirements
final_min_thickness = max(min_thickness_mech, min_thickness_geom)
return final_min_thickness
def generate_grid_pattern(self,
marker_material: Material,
marker_diameter: float,
num_markers: int,
min_edge_distance: float = 5.0) -> List[Marker]:
"""Generates a list of markers arranged in a grid pattern.
Markers are placed within a usable area defined by the minimum edge distance.
The grid attempts to match the panel's aspect ratio.
Args:
marker_material: The material for the markers.
marker_diameter: The diameter for each marker in mm.
num_markers: The total number of markers to generate.
min_edge_distance: Minimum distance from marker center to panel edge in mm.
Returns:
A list of Marker objects positioned in a grid.
"""
if num_markers <= 0:
return []
usable_length = self.length - 2 * min_edge_distance
usable_width = self.width - 2 * min_edge_distance
if usable_length <= 0 or usable_width <= 0:
print("Warning: Panel dimensions too small for minimum edge distance. Placing markers at center.")
center_x = self.length / 2.0
center_y = self.width / 2.0
# Place all markers at the center (they will overlap)
return [Marker(marker_material, marker_diameter, (center_x, center_y)) for _ in range(num_markers)]
# Determine grid dimensions (nx, ny) to be roughly proportional to aspect ratio
aspect_ratio = usable_length / usable_width
nx = max(1, int(np.round(np.sqrt(num_markers * aspect_ratio))))
ny = max(1, int(np.ceil(num_markers / nx)))
# Adjust if grid is too large
while nx * ny < num_markers:
# Increment the smaller dimension ratio relative to aspect ratio
if nx/ny < aspect_ratio:
nx += 1
else:
ny += 1
# Generate grid points within the usable area
x_coords = np.linspace(min_edge_distance, self.length - min_edge_distance, nx)
y_coords = np.linspace(min_edge_distance, self.width - min_edge_distance, ny)
# Create meshgrid and flatten
xv, yv = np.meshgrid(x_coords, y_coords)
grid_points = list(zip(xv.flatten(), yv.flatten()))
# Create markers using the first num_markers grid points
new_markers = []
for i in range(min(num_markers, len(grid_points))):
new_marker = Marker(
material=marker_material,
diameter=marker_diameter,
position=grid_points[i]
)
new_markers.append(new_marker)
if len(new_markers) < num_markers:
print(f"Warning: Could only place {len(new_markers)} markers in the grid. Check parameters.")
return new_markers
def visualize(self, show_stress: bool = False, figsize=(8, 6)):
"""Visualizes the substrate panel and embedded markers.
Args:
show_stress: If True, colors markers based on the simplified stress
concentration factor (Kt). Red indicates high Kt (>=6 or inf),
green indicates low Kt (~1).
figsize: The size of the matplotlib figure.
Returns:
matplotlib figure and axes objects.
"""
fig, ax = plt.subplots(figsize=figsize)
# Draw substrate
ax.add_patch(plt.Rectangle((0, 0), self.length, self.width,
fill=True, alpha=0.2, color='skyblue', label='Substrate'))
# Pre-calculate stress if needed
stress_data = None
if show_stress and self.markers:
stress_data = self.calculate_stress_concentration()
# Draw markers
for i, marker in enumerate(self.markers):
x, y = marker.position
radius = marker.diameter / 2.0
color = 'dimgray' # Default color
label = f'Marker {i}'
if show_stress and stress_data and i in stress_data:
kt = stress_data[i]['stress_concentration_factor']
label += f' (Kt={kt:.2f})'
if kt == float('inf'):
color = 'red'
else:
# Normalize Kt: 1 -> 0 (green), 6+ -> 1 (red)
# Using RdYlGn_r colormap: 0=Red, 0.5=Yellow, 1=Green. We need to invert.
normalized_kt = min(1.0, (kt - 1.0) / 5.0) # Scale Kt=1..6 to 0..1
try:
color = plt.cm.RdYlGn_r(normalized_kt)
except Exception as e:
print(f"Error setting color for Kt={kt}: {e}")
color = 'purple' # Fallback color
circle = Circle((x, y), radius, fill=True, color=color, alpha=0.9)
ax.add_patch(circle)
# Add marker index text inside the circle
ax.text(x, y, str(i), ha='center', va='center', fontsize=max(6, radius*0.8), color='white' if np.mean(plt.cm.colors.to_rgb(color)) < 0.5 else 'black')
# Set axis limits and labels
ax.set_xlim(-0.05 * self.length, 1.05 * self.length)
ax.set_ylim(-0.05 * self.width, 1.05 * self.width)
ax.set_xlabel(f'Length ({self.length} mm)')
ax.set_ylabel(f'Width ({self.width} mm)')
ax.set_title(f'Panel: {self.material.name} ({self.thickness:.1f}mm thick) with {len(self.markers)} markers')
ax.set_aspect('equal', adjustable='box') # Ensure correct aspect ratio
ax.grid(True, linestyle=':', alpha=0.6)
# Add stress colorbar if applicable
if show_stress:
sm = plt.cm.ScalarMappable(cmap=plt.cm.RdYlGn_r, norm=plt.Normalize(vmin=0, vmax=1))
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, label='Simplified Stress Concentration (Kt: 1=Green, >=6=Red)')
cbar.set_ticks([0, 0.25, 0.5, 0.75, 1.0])
cbar.set_ticklabels(['>=6 (High)', '4.75', '3.5', '2.25', '1 (Low)']) # Corresponding Kt values
plt.tight_layout()
return fig, ax
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from typing import List, Dict, Tuple, Optional
import pandas as pd
def setup_matplotlib_style():
"""Set up the matplotlib style for consistent visualizations"""
try:
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette("viridis")
except OSError:
print("Seaborn style 'seaborn-v0_8-whitegrid' not found. Using default.")
plt.style.use('default')
# Optional: set a default palette if seaborn is available
try:
sns.set_palette("viridis")
except NameError:
pass # Seaborn might not be imported if style failed early
def plot_comparison_bars(df: pd.DataFrame,
column: str,
name_column: str,
title: str,
ylabel: str,
ax=None,
rotation: int = 45):
"""Create a bar chart for comparison data
Args:
df: DataFrame containing the data
column: Column name for the y-values
name_column: Column name for the x-values (labels)
title: Plot title
ylabel: Y-axis label
ax: Matplotlib axis to plot on (optional)
rotation: Rotation angle for x-tick labels
Returns:
The matplotlib axis
"""
if ax is None:
_, ax = plt.subplots(figsize=(10, 6))
df.plot(
x=name_column,
y=column,
kind='bar',
ax=ax,
legend=False,
color=sns.color_palette("viridis", len(df))
)
ax.set_title(title)
ax.set_ylabel(ylabel)
ax.tick_params(axis='x', rotation=rotation)
return ax
def plot_substrate_comparisons(substrate_comparison_df: pd.DataFrame):
"""Create comparison visualizations for substrate materials
Args:
substrate_comparison_df: DataFrame with substrate comparison data
"""
if substrate_comparison_df.empty:
print("No data to visualize.")
return
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
plot_comparison_bars(
substrate_comparison_df,
'Radiolucency Score',
'Substrate',
'Radiolucency Comparison',
'Radiolucency Score (higher is better)',
axes[0]
)
plot_comparison_bars(
substrate_comparison_df,
'Min Thickness (mm)',
'Substrate',
'Minimum Thickness Comparison',
'Thickness (mm)',
axes[1]
)
plot_comparison_bars(
substrate_comparison_df,
'Weight (g)',
'Substrate',
'Estimated Weight Comparison',
'Weight (g)',
axes[2]
)
plt.tight_layout()
plt.show()
def plot_marker_comparisons(marker_comparison_df: pd.DataFrame):
"""Create comparison visualizations for marker materials
Args:
marker_comparison_df: DataFrame with marker comparison data
"""
if marker_comparison_df.empty:
print("No data to visualize.")
return
# Bar charts
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
plot_comparison_bars(
marker_comparison_df,
'Radiolucency Score',
'Marker Material',
'Impact on Panel Radiolucency',
'Radiolucency Score (higher is better)',
axes[0]
)
plot_comparison_bars(
marker_comparison_df,
'Relative Visibility Score',
'Marker Material',
'Marker Visibility Comparison',
'Relative Visibility Score (higher is more visible)',
axes[1]
)
plt.tight_layout()
plt.show()
# Scatter plot for trade-off analysis
plt.figure(figsize=(8, 6))
scatter_palette = sns.color_palette("viridis", len(marker_comparison_df))
for i, row in marker_comparison_df.iterrows():
plt.scatter(
row['Radiolucency Score'],
row['Relative Visibility Score'],
s=150,
label=row['Marker Material'],
alpha=0.8,
color=scatter_palette[i]
)
plt.xlabel('Panel Radiolucency Score (higher is better)')
plt.ylabel('Marker Relative Visibility Score (higher is more visible)')
plt.title('Trade-off: Panel Radiolucency vs. Marker Visibility')
plt.grid(True, linestyle=':', alpha=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout(rect=[0, 0, 0.85, 1]) # Adjust layout to make space for legend
plt.show()
print("\nTrade-off Analysis Interpretation:")
print("- Markers in the upper-left quadrant offer the best balance (high visibility, high panel radiolucency).")
print("- Markers in the upper-right have high visibility but significantly reduce panel radiolucency.")
print("- Markers in the lower-left maintain high panel radiolucency but offer low visibility.")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment