|
| 1 | +// MeasureDeviation — directed + symmetric surface deviation (one-sided / |
| 2 | +// symmetric Hausdorff) between two BREPs. |
| 3 | +// |
| 4 | +// Where measure-distance returns the *minimum* gap (≈0 for an overlapping |
| 5 | +// reconstruction-vs-source pair, hence useless as a fidelity figure), this |
| 6 | +// samples one shape's tessellated surface and projects each sample onto the |
| 7 | +// other shape's triangles to report the *worst* / RMS / mean deviation in each |
| 8 | +// direction. The metric a mesh→analytic reconstruction check needs |
| 9 | +// (OCCTMCP #41): fromToTo (over-extension), toToFrom (under-coverage). |
| 10 | +// |
| 11 | +// Mesh-based: both inputs are typically meshes (an STL source vs a STEP |
| 12 | +// reconstruction), and meshing once + a KD-tree is far cheaper than N per-point |
| 13 | +// BRep extrema. The per-sample distance is exact point-to-triangle, so the only |
| 14 | +// approximation is the tessellation — tightened by a finer --deflection. |
| 15 | +// |
| 16 | +// Two input modes: |
| 17 | +// 1. Flag form: |
| 18 | +// occtkit measure-deviation <a.brep> <b.brep> |
| 19 | +// [--deflection D] [--max-samples N] |
| 20 | +// 2. JSON form: |
| 21 | +// { "a": "...", "b": "...", "deflection": D, "maxSamples": N } |
| 22 | + |
| 23 | +import Foundation |
| 24 | +import OCCTSwift |
| 25 | +import ScriptHarness |
| 26 | +import simd |
| 27 | + |
| 28 | +enum MeasureDeviationCommand: Subcommand { |
| 29 | + static let name = "measure-deviation" |
| 30 | + static let summary = "Directed + symmetric surface deviation (Hausdorff) between two BREPs" |
| 31 | + static let usage = """ |
| 32 | + Usage: |
| 33 | + measure-deviation <a.brep> <b.brep> [--deflection D] [--max-samples N] |
| 34 | + measure-deviation <request.json> (JSON request from file) |
| 35 | + measure-deviation (JSON request from stdin) |
| 36 | +
|
| 37 | + --deflection mesh linear deflection (model units). Default: 0.5% of |
| 38 | + the a-shape bbox diagonal. Smaller = finer = tighter bound. |
| 39 | + --max-samples max source surface samples per direction. Default 20000. |
| 40 | + """ |
| 41 | + |
| 42 | + private struct Request { |
| 43 | + var a: String |
| 44 | + var b: String |
| 45 | + var deflection: Double? |
| 46 | + var maxSamples: Int |
| 47 | + } |
| 48 | + |
| 49 | + private struct JSONRequest: Decodable { |
| 50 | + let a: String |
| 51 | + let b: String |
| 52 | + let deflection: Double? |
| 53 | + let maxSamples: Int? |
| 54 | + } |
| 55 | + |
| 56 | + struct DirectionStat: Encodable { |
| 57 | + let max: Double |
| 58 | + let rms: Double |
| 59 | + let mean: Double |
| 60 | + let worstPoint: [Double] |
| 61 | + let samples: Int |
| 62 | + } |
| 63 | + |
| 64 | + struct Response: Encodable { |
| 65 | + let deflection: Double |
| 66 | + let fromToTo: DirectionStat |
| 67 | + let toToFrom: DirectionStat |
| 68 | + let symmetricHausdorff: Double |
| 69 | + } |
| 70 | + |
| 71 | + static func run(args: [String]) throws -> Int32 { |
| 72 | + let req = try parseRequest(args: args) |
| 73 | + let aShape = try GraphIO.loadBREP(at: req.a) |
| 74 | + let bShape = try GraphIO.loadBREP(at: req.b) |
| 75 | + |
| 76 | + let defl = req.deflection ?? defaultDeflection(for: aShape) |
| 77 | + guard defl > 0 else { throw ScriptError.message("deflection must be positive") } |
| 78 | + |
| 79 | + guard let aTris = TriMesh(shape: aShape, deflection: defl) else { |
| 80 | + throw ScriptError.message("Failed to tessellate '\(req.a)' for deviation") |
| 81 | + } |
| 82 | + guard let bTris = TriMesh(shape: bShape, deflection: defl) else { |
| 83 | + throw ScriptError.message("Failed to tessellate '\(req.b)' for deviation") |
| 84 | + } |
| 85 | + guard let fwd = directedDeviation(source: aTris, target: bTris, maxSamples: req.maxSamples), |
| 86 | + let rev = directedDeviation(source: bTris, target: aTris, maxSamples: req.maxSamples) else { |
| 87 | + throw ScriptError.message("Deviation computation failed (empty tessellation)") |
| 88 | + } |
| 89 | + |
| 90 | + try GraphIO.emitJSON(Response( |
| 91 | + deflection: defl, |
| 92 | + fromToTo: fwd, |
| 93 | + toToFrom: rev, |
| 94 | + symmetricHausdorff: Swift.max(fwd.max, rev.max) |
| 95 | + )) |
| 96 | + return 0 |
| 97 | + } |
| 98 | + |
| 99 | + // ── tessellation snapshot ─────────────────────────────────────────── |
| 100 | + |
| 101 | + struct TriMesh { |
| 102 | + let vertices: [SIMD3<Double>] |
| 103 | + let triangles: [(UInt32, UInt32, UInt32)] |
| 104 | + let kd: KDTree |
| 105 | + let incident: [[Int]] |
| 106 | + |
| 107 | + init?(shape: Shape, deflection: Double) { |
| 108 | + var params = MeshParameters.default |
| 109 | + params.deflection = deflection |
| 110 | + params.internalVertices = true |
| 111 | + params.inParallel = true |
| 112 | + params.allowQualityDecrease = true // honour the requested deflection (#211) |
| 113 | + guard let mesh = shape.mesh(parameters: params) else { return nil } |
| 114 | + |
| 115 | + let verts = mesh.vertices.map { SIMD3<Double>(Double($0.x), Double($0.y), Double($0.z)) } |
| 116 | + let idx = mesh.indices |
| 117 | + guard !verts.isEmpty, idx.count >= 3, let kd = KDTree(points: verts) else { return nil } |
| 118 | + |
| 119 | + var tris: [(UInt32, UInt32, UInt32)] = [] |
| 120 | + tris.reserveCapacity(idx.count / 3) |
| 121 | + var adj = [[Int]](repeating: [], count: verts.count) |
| 122 | + var t = 0 |
| 123 | + while t + 2 < idx.count { |
| 124 | + let a = idx[t], b = idx[t + 1], c = idx[t + 2] |
| 125 | + let ti = tris.count |
| 126 | + tris.append((a, b, c)) |
| 127 | + adj[Int(a)].append(ti); adj[Int(b)].append(ti); adj[Int(c)].append(ti) |
| 128 | + t += 3 |
| 129 | + } |
| 130 | + self.vertices = verts; self.triangles = tris; self.kd = kd; self.incident = adj |
| 131 | + } |
| 132 | + } |
| 133 | + |
| 134 | + // ── directed deviation ────────────────────────────────────────────── |
| 135 | + |
| 136 | + static func directedDeviation(source: TriMesh, target: TriMesh, maxSamples: Int) -> DirectionStat? { |
| 137 | + let n = source.vertices.count |
| 138 | + guard n > 0 else { return nil } |
| 139 | + let stride = maxSamples > 0 ? Swift.max(1, (n + maxSamples - 1) / maxSamples) : 1 |
| 140 | + let k = 6 |
| 141 | + |
| 142 | + var maxD = 0.0, sumSq = 0.0, sum = 0.0, count = 0 |
| 143 | + var worst = SIMD3<Double>(0, 0, 0) |
| 144 | + var stamp = [Int](repeating: -1, count: target.triangles.count) |
| 145 | + |
| 146 | + var i = 0 |
| 147 | + while i < n { |
| 148 | + let p = source.vertices[i] |
| 149 | + var best = Double.greatestFiniteMagnitude |
| 150 | + let neighbours = target.kd.kNearest(to: p, k: k) |
| 151 | + if neighbours.isEmpty { |
| 152 | + if let nv = target.kd.nearest(to: p) { best = nv.distance } |
| 153 | + } else { |
| 154 | + for (vi, _) in neighbours { |
| 155 | + for ti in target.incident[vi] where stamp[ti] != i { |
| 156 | + stamp[ti] = i |
| 157 | + let (a, b, c) = target.triangles[ti] |
| 158 | + let d = pointTriangleDistance(p, target.vertices[Int(a)], |
| 159 | + target.vertices[Int(b)], target.vertices[Int(c)]) |
| 160 | + if d < best { best = d } |
| 161 | + } |
| 162 | + } |
| 163 | + if best == .greatestFiniteMagnitude, let nv = target.kd.nearest(to: p) { best = nv.distance } |
| 164 | + } |
| 165 | + if best != .greatestFiniteMagnitude { |
| 166 | + if best > maxD { maxD = best; worst = p } |
| 167 | + sumSq += best * best; sum += best; count += 1 |
| 168 | + } |
| 169 | + i += stride |
| 170 | + } |
| 171 | + guard count > 0 else { return nil } |
| 172 | + return DirectionStat( |
| 173 | + max: maxD, |
| 174 | + rms: (sumSq / Double(count)).squareRoot(), |
| 175 | + mean: sum / Double(count), |
| 176 | + worstPoint: [worst.x, worst.y, worst.z], |
| 177 | + samples: count |
| 178 | + ) |
| 179 | + } |
| 180 | + |
| 181 | + // ── geometry helpers ──────────────────────────────────────────────── |
| 182 | + |
| 183 | + static func defaultDeflection(for shape: Shape) -> Double { |
| 184 | + let b = shape.bounds |
| 185 | + let diag = simd_length(b.max - b.min) |
| 186 | + return Swift.max(diag * 0.005, 1e-6) |
| 187 | + } |
| 188 | + |
| 189 | + /// Closest-point-on-triangle distance, Ericson "Real-Time Collision |
| 190 | + /// Detection" §5.1.5. |
| 191 | + static func pointTriangleDistance(_ p: SIMD3<Double>, _ a: SIMD3<Double>, |
| 192 | + _ b: SIMD3<Double>, _ c: SIMD3<Double>) -> Double { |
| 193 | + let ab = b - a, ac = c - a, ap = p - a |
| 194 | + let d1 = simd_dot(ab, ap), d2 = simd_dot(ac, ap) |
| 195 | + if d1 <= 0 && d2 <= 0 { return simd_length(ap) } |
| 196 | + let bp = p - b |
| 197 | + let d3 = simd_dot(ab, bp), d4 = simd_dot(ac, bp) |
| 198 | + if d3 >= 0 && d4 <= d3 { return simd_length(bp) } |
| 199 | + let vc = d1 * d4 - d3 * d2 |
| 200 | + if vc <= 0 && d1 >= 0 && d3 <= 0 { let v = d1 / (d1 - d3); return simd_length(p - (a + v * ab)) } |
| 201 | + let cp = p - c |
| 202 | + let d5 = simd_dot(ab, cp), d6 = simd_dot(ac, cp) |
| 203 | + if d6 >= 0 && d5 <= d6 { return simd_length(cp) } |
| 204 | + let vb = d5 * d2 - d1 * d6 |
| 205 | + if vb <= 0 && d2 >= 0 && d6 <= 0 { let w = d2 / (d2 - d6); return simd_length(p - (a + w * ac)) } |
| 206 | + let va = d3 * d6 - d5 * d4 |
| 207 | + if va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0 { |
| 208 | + let w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); return simd_length(p - (b + w * (c - b))) |
| 209 | + } |
| 210 | + let denom = 1.0 / (va + vb + vc) |
| 211 | + return simd_length(p - (a + ab * (vb * denom) + ac * (vc * denom))) |
| 212 | + } |
| 213 | + |
| 214 | + // ── request parsing ───────────────────────────────────────────────── |
| 215 | + |
| 216 | + private static func parseRequest(args: [String]) throws -> Request { |
| 217 | + if args.count == 1, args[0].hasSuffix(".json") { |
| 218 | + return try decodeJSON(data: try readFile(args[0])) |
| 219 | + } |
| 220 | + if args.isEmpty { |
| 221 | + return try decodeJSON(data: FileHandle.standardInput.readDataToEndOfFile()) |
| 222 | + } |
| 223 | + guard args.count >= 2, !args[0].hasPrefix("-"), !args[1].hasPrefix("-") else { |
| 224 | + throw ScriptError.message("Expected: <a.brep> <b.brep> [flags]") |
| 225 | + } |
| 226 | + var deflection: Double? |
| 227 | + var maxSamples = 20_000 |
| 228 | + var i = 2 |
| 229 | + while i < args.count { |
| 230 | + switch args[i] { |
| 231 | + case "--deflection": |
| 232 | + i += 1 |
| 233 | + guard i < args.count, let d = Double(args[i]) else { |
| 234 | + throw ScriptError.message("--deflection expects a number") |
| 235 | + } |
| 236 | + deflection = d |
| 237 | + case "--max-samples": |
| 238 | + i += 1 |
| 239 | + guard i < args.count, let n = Int(args[i]) else { |
| 240 | + throw ScriptError.message("--max-samples expects an integer") |
| 241 | + } |
| 242 | + maxSamples = n |
| 243 | + default: |
| 244 | + throw ScriptError.message("Unknown flag: \(args[i])") |
| 245 | + } |
| 246 | + i += 1 |
| 247 | + } |
| 248 | + return Request(a: args[0], b: args[1], deflection: deflection, maxSamples: maxSamples) |
| 249 | + } |
| 250 | + |
| 251 | + private static func readFile(_ path: String) throws -> Data { |
| 252 | + guard let bytes = FileManager.default.contents(atPath: path) else { |
| 253 | + throw ScriptError.message("Failed to read request at \(path)") |
| 254 | + } |
| 255 | + return bytes |
| 256 | + } |
| 257 | + |
| 258 | + private static func decodeJSON(data: Data) throws -> Request { |
| 259 | + let raw: JSONRequest |
| 260 | + do { |
| 261 | + raw = try JSONDecoder().decode(JSONRequest.self, from: data) |
| 262 | + } catch { |
| 263 | + throw ScriptError.message("Invalid JSON: \(error.localizedDescription)") |
| 264 | + } |
| 265 | + return Request(a: raw.a, b: raw.b, deflection: raw.deflection, |
| 266 | + maxSamples: raw.maxSamples ?? 20_000) |
| 267 | + } |
| 268 | +} |
0 commit comments