add linear/circle Nyquist fit resolver to iOS
This commit is contained in:
parent
6a09782d30
commit
1abf46f0c3
|
|
@ -285,46 +285,68 @@ struct EISView: View {
|
|||
with: .color(nyqColor))
|
||||
}
|
||||
|
||||
// circle fit
|
||||
if let fit = kasaCircleFit(points: points.map { (Double($0.zReal), Double(-$0.zImag)) }) {
|
||||
let disc = fit.r * fit.r - fit.cy * fit.cy
|
||||
if disc > 0 {
|
||||
let sd = sqrt(disc)
|
||||
let rs = fit.cx - sd
|
||||
let rp = 2 * sd
|
||||
// fit overlay
|
||||
let fitColor = Color(red: 0.7, green: 0.3, blue: 0.9).opacity(0.5)
|
||||
let fitPtColor = Color(red: 0.9, green: 0.4, blue: 0.9)
|
||||
|
||||
if rp > 0 {
|
||||
let thetaR = atan2(-fit.cy, sd)
|
||||
var thetaL = atan2(-fit.cy, -sd)
|
||||
if thetaL < thetaR { thetaL += 2 * .pi }
|
||||
if let result = fitNyquist(points: points.map { (Double($0.zReal), Double(-$0.zImag)) }) {
|
||||
switch result {
|
||||
case .circle(let fit):
|
||||
let disc = fit.r * fit.r - fit.cy * fit.cy
|
||||
if disc > 0 {
|
||||
let sd = sqrt(disc)
|
||||
let rs = fit.cx - sd
|
||||
let rp = 2 * sd
|
||||
|
||||
let nArc = 120
|
||||
var arcPath = Path()
|
||||
for i in 0...nArc {
|
||||
let t = thetaR + (thetaL - thetaR) * Double(i) / Double(nArc)
|
||||
let ax = fit.cx + fit.r * cos(t)
|
||||
let ay = fit.cy + fit.r * sin(t)
|
||||
let pt = CGPoint(x: lx(CGFloat(ax)), y: ly(CGFloat(ay)))
|
||||
if i == 0 { arcPath.move(to: pt) } else { arcPath.addLine(to: pt) }
|
||||
if rp > 0 {
|
||||
let thetaR = atan2(-fit.cy, sd)
|
||||
var thetaL = atan2(-fit.cy, -sd)
|
||||
if thetaL < thetaR { thetaL += 2 * .pi }
|
||||
|
||||
let nArc = 120
|
||||
var arcPath = Path()
|
||||
for i in 0...nArc {
|
||||
let t = thetaR + (thetaL - thetaR) * Double(i) / Double(nArc)
|
||||
let ax = fit.cx + fit.r * cos(t)
|
||||
let ay = fit.cy + fit.r * sin(t)
|
||||
let pt = CGPoint(x: lx(CGFloat(ax)), y: ly(CGFloat(ay)))
|
||||
if i == 0 { arcPath.move(to: pt) } else { arcPath.addLine(to: pt) }
|
||||
}
|
||||
context.stroke(arcPath, with: .color(fitColor), lineWidth: 1.5)
|
||||
|
||||
let rsScr = CGPoint(x: lx(CGFloat(rs)), y: ly(0))
|
||||
let rpScr = CGPoint(x: lx(CGFloat(rs + rp)), y: ly(0))
|
||||
context.fill(Path(ellipseIn: CGRect(x: rsScr.x - 5, y: rsScr.y - 5, width: 10, height: 10)),
|
||||
with: .color(fitPtColor))
|
||||
context.fill(Path(ellipseIn: CGRect(x: rpScr.x - 5, y: rpScr.y - 5, width: 10, height: 10)),
|
||||
with: .color(fitPtColor))
|
||||
context.draw(
|
||||
Text(String(format: "Rs=%.0f", rs)).font(.caption2).foregroundStyle(fitPtColor),
|
||||
at: CGPoint(x: rsScr.x, y: rsScr.y + 14))
|
||||
context.draw(
|
||||
Text(String(format: "Rp=%.0f", rp)).font(.caption2).foregroundStyle(fitPtColor),
|
||||
at: CGPoint(x: rpScr.x, y: rpScr.y + 14))
|
||||
}
|
||||
context.stroke(arcPath, with: .color(Color(red: 0.7, green: 0.3, blue: 0.9).opacity(0.5)),
|
||||
lineWidth: 1.5)
|
||||
|
||||
// Rs and Rp markers
|
||||
let fitPtColor = Color(red: 0.9, green: 0.4, blue: 0.9)
|
||||
let rsScr = CGPoint(x: lx(CGFloat(rs)), y: ly(0))
|
||||
let rpScr = CGPoint(x: lx(CGFloat(rs + rp)), y: ly(0))
|
||||
context.fill(Path(ellipseIn: CGRect(x: rsScr.x - 5, y: rsScr.y - 5, width: 10, height: 10)),
|
||||
with: .color(fitPtColor))
|
||||
context.fill(Path(ellipseIn: CGRect(x: rpScr.x - 5, y: rpScr.y - 5, width: 10, height: 10)),
|
||||
with: .color(fitPtColor))
|
||||
context.draw(
|
||||
Text(String(format: "Rs=%.0f", rs)).font(.caption2).foregroundStyle(fitPtColor),
|
||||
at: CGPoint(x: rsScr.x, y: rsScr.y + 14))
|
||||
context.draw(
|
||||
Text(String(format: "Rp=%.0f", rp)).font(.caption2).foregroundStyle(fitPtColor),
|
||||
at: CGPoint(x: rpScr.x, y: rpScr.y + 14))
|
||||
}
|
||||
case .linear(let fit):
|
||||
let xVals = points.map { CGFloat($0.zReal) }.filter { $0.isFinite }
|
||||
guard let xMin = xVals.min(), let xMax = xVals.max() else { break }
|
||||
let pad = (xMax - xMin) * 0.1
|
||||
let x0 = Double(xMin - pad)
|
||||
let x1 = Double(xMax + pad)
|
||||
let y0 = fit.slope * x0 + fit.yIntercept
|
||||
let y1 = fit.slope * x1 + fit.yIntercept
|
||||
var linePath = Path()
|
||||
linePath.move(to: CGPoint(x: lx(CGFloat(x0)), y: ly(CGFloat(y0))))
|
||||
linePath.addLine(to: CGPoint(x: lx(CGFloat(x1)), y: ly(CGFloat(y1))))
|
||||
context.stroke(linePath, with: .color(fitColor), lineWidth: 1.5)
|
||||
|
||||
let rsScr = CGPoint(x: lx(CGFloat(fit.rs)), y: ly(0))
|
||||
context.fill(Path(ellipseIn: CGRect(x: rsScr.x - 5, y: rsScr.y - 5, width: 10, height: 10)),
|
||||
with: .color(fitPtColor))
|
||||
context.draw(
|
||||
Text(String(format: "Rs=%.0f", fit.rs)).font(.caption2).foregroundStyle(fitPtColor),
|
||||
at: CGPoint(x: rsScr.x, y: rsScr.y + 14))
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -370,7 +392,7 @@ struct EISView: View {
|
|||
}
|
||||
}
|
||||
|
||||
// MARK: - Kasa circle fit (ported from plot.rs)
|
||||
// MARK: - Nyquist fit
|
||||
|
||||
struct CircleFitResult {
|
||||
let cx: Double
|
||||
|
|
@ -378,6 +400,17 @@ struct CircleFitResult {
|
|||
let r: Double
|
||||
}
|
||||
|
||||
struct LinearFitResult {
|
||||
let slope: Double
|
||||
let yIntercept: Double
|
||||
let rs: Double
|
||||
}
|
||||
|
||||
enum NyquistFitResult {
|
||||
case circle(CircleFitResult)
|
||||
case linear(LinearFitResult)
|
||||
}
|
||||
|
||||
func kasaCircleFit(points: [(Double, Double)]) -> CircleFitResult? {
|
||||
let all = points.filter { $0.0.isFinite && $0.1.isFinite }
|
||||
guard all.count >= 4 else { return nil }
|
||||
|
|
@ -463,6 +496,57 @@ private func kasaFitRaw(_ pts: [(Double, Double)]) -> (Double, Double, Double)?
|
|||
return (cx, cy, sqrt(rSq))
|
||||
}
|
||||
|
||||
private func cumulativeTurning(_ pts: [(Double, Double)]) -> Double {
|
||||
guard pts.count >= 3 else { return 0 }
|
||||
var total = 0.0
|
||||
for i in 1..<(pts.count - 1) {
|
||||
let (dx1, dy1) = (pts[i].0 - pts[i-1].0, pts[i].1 - pts[i-1].1)
|
||||
let (dx2, dy2) = (pts[i+1].0 - pts[i].0, pts[i+1].1 - pts[i].1)
|
||||
let cross = dx1 * dy2 - dy1 * dx2
|
||||
let dot = dx1 * dx2 + dy1 * dy2
|
||||
total += abs(atan2(cross, dot))
|
||||
}
|
||||
return total
|
||||
}
|
||||
|
||||
private func fitLinear(_ pts: [(Double, Double)]) -> LinearFitResult? {
|
||||
guard pts.count >= 2 else { return nil }
|
||||
let n = Double(pts.count)
|
||||
let sx = pts.map(\.0).reduce(0, +)
|
||||
let sy = pts.map(\.1).reduce(0, +)
|
||||
let sx2 = pts.map { $0.0 * $0.0 }.reduce(0, +)
|
||||
let sxy = pts.map { $0.0 * $0.1 }.reduce(0, +)
|
||||
let denom = n * sx2 - sx * sx
|
||||
guard abs(denom) > 1e-20 else { return nil }
|
||||
let slope = (n * sxy - sx * sy) / denom
|
||||
let yInt = (sy - slope * sx) / n
|
||||
let rs = abs(slope) > 1e-10 ? -yInt / slope : sx / n
|
||||
return LinearFitResult(slope: slope, yIntercept: yInt, rs: rs)
|
||||
}
|
||||
|
||||
func fitNyquist(points: [(Double, Double)]) -> NyquistFitResult? {
|
||||
let all = points.filter { $0.0.isFinite && $0.1.isFinite }
|
||||
guard all.count >= 4 else { return nil }
|
||||
|
||||
if cumulativeTurning(all) < 0.524 {
|
||||
if let lin = fitLinear(all) { return .linear(lin) }
|
||||
}
|
||||
|
||||
if let circle = kasaCircleFit(points: points) {
|
||||
let avgErr = all.map { p in
|
||||
abs(sqrt((p.0 - circle.cx) * (p.0 - circle.cx) +
|
||||
(p.1 - circle.cy) * (p.1 - circle.cy)) - circle.r)
|
||||
}.reduce(0, +) / (Double(all.count) * circle.r)
|
||||
if avgErr > 0.15 {
|
||||
if let lin = fitLinear(all) { return .linear(lin) }
|
||||
}
|
||||
return .circle(circle)
|
||||
}
|
||||
|
||||
if let lin = fitLinear(all) { return .linear(lin) }
|
||||
return nil
|
||||
}
|
||||
|
||||
// MARK: - Canvas drawing helpers
|
||||
|
||||
private func drawPolyline(context: GraphicsContext, points: [CGPoint], color: Color, width: CGFloat) {
|
||||
|
|
|
|||
Loading…
Reference in New Issue