//! Ablation zone computation and isosurface extraction //! //! Computes the ablation zone (tissue above 60°C) and generates //! mesh data for visualization. use bioheat_shared::{AblationZone, TemperatureField}; use serde::{Deserialize, Serialize}; /// Ablation threshold temperature in °C pub const ABLATION_THRESHOLD: f32 = 60.0; /// Computes ablation zone from temperature field #[derive(Debug, Clone, Serialize, Deserialize)] pub struct AblationZoneComputer { /// Temperature threshold for ablation (typically 60°C) pub threshold: f32, } impl AblationZoneComputer { /// Create new computer with default threshold #[must_use] pub fn new() -> Self { Self { threshold: ABLATION_THRESHOLD, } } /// Create with custom threshold #[must_use] pub fn with_threshold(threshold: f32) -> Self { Self { threshold } } /// Compute ablation zone from temperature field #[must_use] pub fn compute(&self, field: &TemperatureField) -> AblationZone { let volume = self.compute_volume(field); let dimensions = self.compute_dimensions(field); // For now, we don't generate the full isosurface mesh // (marching cubes is complex and would add significant code) // We just compute volume and dimensions AblationZone { vertices: Vec::new(), indices: Vec::new(), volume, volume_mm3: volume * 1e9, // m³ to mm³ threshold: self.threshold, dimensions, } } /// Compute ablated volume in m³ fn compute_volume(&self, field: &TemperatureField) -> f32 { let spacing = field.spacing(); let voxel_volume = spacing.x * spacing.y * spacing.z; let ablated_voxels = field .values .iter() .filter(|&&t| t >= self.threshold) .count(); ablated_voxels as f32 * voxel_volume } /// Compute approximate dimensions of ablation zone (bounding box) fn compute_dimensions(&self, field: &TemperatureField) -> (f32, f32, f32) { let (nx, ny, nz) = field.resolution; let spacing = field.spacing(); // Find extent in each dimension let mut x_min = nx; let mut x_max = 0; let mut y_min = ny; let mut y_max = 0; let mut z_min = nz; let mut z_max = 0; for k in 0..nz { for j in 0..ny { for i in 0..nx { if field.at(i, j, k) >= self.threshold { x_min = x_min.min(i); x_max = x_max.max(i); y_min = y_min.min(j); y_max = y_max.max(j); z_min = z_min.min(k); z_max = z_max.max(k); } } } } // Convert to physical dimensions if x_max >= x_min { ( (x_max - x_min + 1) as f32 * spacing.x, (y_max - y_min + 1) as f32 * spacing.y, (z_max - z_min + 1) as f32 * spacing.z, ) } else { (0.0, 0.0, 0.0) // No ablation } } /// Check if ablation zone meets minimum size requirements #[must_use] pub fn meets_minimum_size(&self, zone: &AblationZone, min_volume_mm3: f32) -> bool { zone.volume_mm3 >= min_volume_mm3 } /// Estimate the equivalent sphere diameter of ablation zone #[must_use] pub fn equivalent_sphere_diameter_mm(volume_mm3: f32) -> f32 { // V = (4/3)πr³ => d = 2 * (3V/(4π))^(1/3) let r_mm = (3.0 * volume_mm3 / (4.0 * std::f32::consts::PI)).powf(1.0 / 3.0); 2.0 * r_mm } } impl Default for AblationZoneComputer { fn default() -> Self { Self::new() } } /// Statistics about the ablation zone #[derive(Debug, Clone, Serialize, Deserialize)] pub struct AblationStatistics { /// Volume in mm³ pub volume_mm3: f32, /// Equivalent sphere diameter in mm pub equivalent_diameter_mm: f32, /// Dimensions (width, height, depth) in mm pub dimensions_mm: (f32, f32, f32), /// Maximum temperature in the ablation zone pub max_temperature: f32, /// Percentage of domain ablated pub ablation_percentage: f32, } impl AblationStatistics { /// Compute statistics from temperature field #[must_use] pub fn from_field(field: &TemperatureField, threshold: f32) -> Self { let computer = AblationZoneComputer::with_threshold(threshold); let zone = computer.compute(field); let domain_volume = field.bounds.volume() * 1e9; // m³ to mm³ let ablation_percentage = if domain_volume > 0.0 { 100.0 * zone.volume_mm3 / domain_volume } else { 0.0 }; Self { volume_mm3: zone.volume_mm3, equivalent_diameter_mm: AblationZoneComputer::equivalent_sphere_diameter_mm( zone.volume_mm3, ), dimensions_mm: ( zone.dimensions.0 * 1000.0, zone.dimensions.1 * 1000.0, zone.dimensions.2 * 1000.0, ), max_temperature: field.max_temperature(), ablation_percentage, } } } #[cfg(test)] mod tests { use super::*; use bioheat_shared::BoundingBox3D; fn create_test_field(with_hot_spot: bool) -> TemperatureField { let bounds = BoundingBox3D::from_dimensions(0.1, 0.1, 0.1); let resolution = (10, 10, 10); let mut field = TemperatureField::body_temperature(resolution, bounds); if with_hot_spot { // Create a hot region in the center for k in 4..7 { for j in 4..7 { for i in 4..7 { field.set(i, j, k, 80.0); // Above ablation threshold } } } } field } #[test] fn test_no_ablation() { let field = create_test_field(false); let computer = AblationZoneComputer::new(); let zone = computer.compute(&field); assert!(zone.volume < 1e-12); assert!(zone.is_empty()); } #[test] fn test_with_ablation() { let field = create_test_field(true); let computer = AblationZoneComputer::new(); let zone = computer.compute(&field); assert!(zone.volume > 0.0); assert!(zone.volume_mm3 > 0.0); assert!(!zone.is_empty()); } #[test] fn test_dimensions() { let field = create_test_field(true); let computer = AblationZoneComputer::new(); let zone = computer.compute(&field); let (w, h, d) = zone.dimensions; assert!(w > 0.0); assert!(h > 0.0); assert!(d > 0.0); } #[test] fn test_equivalent_diameter() { // 1 mL = 1000 mm³ should give ~12.4 mm diameter let diameter = AblationZoneComputer::equivalent_sphere_diameter_mm(1000.0); assert!((diameter - 12.4).abs() < 0.2); } #[test] fn test_statistics() { let field = create_test_field(true); let stats = AblationStatistics::from_field(&field, ABLATION_THRESHOLD); assert!(stats.volume_mm3 > 0.0); assert!(stats.max_temperature > 60.0); assert!(stats.ablation_percentage > 0.0); assert!(stats.ablation_percentage < 100.0); } }