refactor: BV estimation — y-z/x-z ellipse fitting, wall repair, outlier removal

Port of appshare #21 merge (bv_estimation.py):
- repairCenterWallsFor6ch: gap=1 보간, gap>=2 top 그룹 제거
- computeLrRatio: chord 비율 → x-z 6점 타원 피팅 (b_lr/a_ap)
- Post-median outlier 제거 (median ±25%, ≥3ch)
- Cap 높이: ellipseCapHeights 삭제 → y-z 8점 타원 피팅 inline
- Top cap 상한: 미검출 상위 채널 빔 y좌표 기준 제한

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

app/src/main/java/com/example/medilightv2android/managers/PiezoBVEstimator.kt

@@ -112,33 +112,69 @@ object PiezoHW {
         }
 }
 
-// ── LR Ratio Computation (from lateral CH4/CH5) ──
+// ── Center Wall Repair (6ch) ──
 
+/**
+ * 6ch center 채널(CH0~CH3) 패턴을 BV 계산용으로 보정.
+ * - gap 1개: 선형 보간
+ * - gap 2개+: 위쪽 그룹 버리고 아래쪽만 사용
+ * Port of bv_estimation.py _repair_center_walls_for_6ch (#21 merge)
+ */
+private fun repairCenterWallsFor6ch(
+    centerWalls: List<Pair<Int, Int>?>
+): List<Pair<Int, Int>?> {
+    val result = centerWalls.toMutableList()
+    val valid = result.indices.filter { i ->
+        val w = result[i]; w != null
+    }
+    if (valid.size < 2) return result
+
+    for (k in 0 until valid.size - 1) {
+        val prevIdx = valid[k]
+        val nextIdx = valid[k + 1]
+        val gap = nextIdx - prevIdx - 1
+        if (gap <= 0) continue
+        if (gap == 1) {
+            val pw = result[prevIdx]!!
+            val nw = result[nextIdx]!!
+            val ant = ((pw.first + nw.first) / 2.0 + 0.5).toInt()
+            val post = ((pw.second + nw.second) / 2.0 + 0.5).toInt()
+            result[prevIdx + 1] = Pair(ant, post)
+            continue
+        }
+        // gap >= 2: drop top group
+        for (i in 0 until nextIdx) result[i] = null
+        return result
+    }
+    return result
+}
+
+// ── LR Ratio Computation (x-z ellipse fitting, #21 merge) ──
+
+/**
+ * x-z 평면 타원 피팅 기반 LR/AP ratio.
+ * CH4/CH5의 SI 높이에서 center(CH1,CH2) 벽 좌표를 보간하고,
+ * lateral 벽 좌표와 합쳐 6개 경계점으로 x-z 평면 타원 피팅.
+ * lr_ratio = LR 반축 / AP 반축.
+ */
 fun computeLrRatio(
     centerWalls: List<Pair<Int, Int>?>,
     lateralWalls: List<Pair<Int, Int>?>,
-    maxRatio: Double = 1.0
+    maxRatio: Double = 1.0,
+    sensorZMm: DoubleArray = PiezoHW.sensorZMmAll
 ): Double {
     val hw = PiezoHW
     val neighbors = hw.lateralNeighbors
-
-    // Center AP chord from neighbors
-    val centerDs = mutableListOf<Double>()
-    for (ni in neighbors) {
-        if (ni < centerWalls.size) {
-            val w = centerWalls[ni] ?: continue
-            val (dNear, dFar) = segmentToDistancesMm(w.first, w.second)
-            val theta = hw.degreeAll[hw.centerCh[ni]] * Math.PI / 180.0
-            centerDs.add((dFar - dNear) * cos(theta))
-        }
+    val availableNbrs = neighbors.filter { ni ->
+        ni < centerWalls.size && centerWalls[ni] != null
     }
-    if (centerDs.isEmpty()) return hw.lrRatioNoDetection
-    val dCenter = centerDs.average()
-    if (dCenter <= 0) return hw.lrRatioInvalid
+    if (availableNbrs.size < 2) return hw.lrRatioNoDetection
 
-    // Lateral chords
-    data class LatResult(val ratio: Double, val yMid: Double)
-    val latResults = mutableListOf<LatResult>()
+    // Lateral wall points (x-z plane)
+    val latPtsX = mutableListOf<Double>()
+    val latPtsZ = mutableListOf<Double>()
+    val latSiPositions = mutableListOf<Double>()
+    val latChords = mutableListOf<Double>()
 
     for ((li, lCh) in hw.lateralCh.withIndex()) {
         if (li >= lateralWalls.size) continue
@@ -147,59 +183,72 @@ fun computeLrRatio(
         val alpha = hw.degreeAll[lCh] * Math.PI / 180.0
         val beta = hw.degreeLRAll[lCh] * Math.PI / 180.0
         val sx = hw.sensorXMmAll[lCh]
+        for (d in listOf(dNear, dFar)) {
+            latPtsX.add(sx + d * cos(alpha) * sin(beta))
+            latPtsZ.add(d * cos(alpha) * cos(beta))
+        }
+        latSiPositions.add(sensorZMm[lCh])
+        latChords.add(abs(dFar - dNear))
+    }
+    if (latPtsX.isEmpty()) return hw.lrRatioNoDetection
+    val targetSi = latSiPositions.average()
 
-        val pFactor = cos(alpha) * cos(beta)
-        val qFactor = cos(alpha) * sin(beta)
-
-        val dLateral = (dFar - dNear) * pFactor
-        val dMid = (dNear + dFar) / 2.0
-        val yMid = sx + dMid * qFactor
-
-        val ratio = dLateral / dCenter
-        if (ratio <= 0 || ratio >= maxRatio) continue
-        if (abs(yMid) < 1e-3) continue
-        latResults.add(LatResult(ratio, yMid))
+    // Center wall z-coords at neighbor channels, interpolated to lateral SI height
+    val nbrSi = mutableListOf<Double>()
+    val nbrZAnt = mutableListOf<Double>()
+    val nbrZPost = mutableListOf<Double>()
+    for (ni in availableNbrs) {
+        val w = centerWalls[ni]!!
+        val ch = hw.centerCh[ni]
+        val (dNear, dFar) = segmentToDistancesMm(w.first, w.second)
+        val th = hw.degreeAll[ch] * Math.PI / 180.0
+        nbrSi.add(sensorZMm[ch])
+        nbrZAnt.add(dNear * cos(th))
+        nbrZPost.add(dFar * cos(th))
     }
 
-    if (latResults.isEmpty()) return hw.lrRatioNoDetection
-
-    val lrRaw: Double
-    if (latResults.size >= 2) {
-        // Two laterals: solve ellipse
-        val yL = latResults[0].yMid
-        val rL = latResults[0].ratio
-        val yR = latResults[1].yMid
-        val rR = latResults[1].ratio
-
-        val denom = yR * (rL * rL - 1.0) - yL * (rR * rR - 1.0)
-        if (abs(denom) < 1e-12) return hw.lrRatioInvalid
-
-        val u = yL * yR * (yR - yL) / denom
-        if (u <= 0) return hw.lrRatioInvalid
-
-        val d = -(u * (rL * rL - 1.0) + yL * yL) / (2.0 * yL)
-        val bSq = u + d * d
-        if (bSq <= 0) return hw.lrRatioInvalid
-
-        lrRaw = 2.0 * sqrt(u) / dCenter
+    val zAntCenter: Double
+    val zPostCenter: Double
+    if (nbrSi.size >= 2) {
+        val w = if (abs(nbrSi[0] - nbrSi[1]) > 1e-6)
+            ((targetSi - nbrSi[1]) / (nbrSi[0] - nbrSi[1])).coerceIn(0.0, 1.0) else 0.5
+        zAntCenter = nbrZAnt[1] + w * (nbrZAnt[0] - nbrZAnt[1])
+        zPostCenter = nbrZPost[1] + w * (nbrZPost[0] - nbrZPost[1])
     } else {
-        // Single lateral: assume center
-        val r = latResults[0].ratio
-        val yMid = latResults[0].yMid
-        if (r >= 1.0) return hw.lrRatioInvalid
-        val bEst = abs(yMid) / sqrt(1.0 - r * r)
-        val aEst = dCenter / 2.0
-        lrRaw = bEst / aEst
+        zAntCenter = nbrZAnt[0]
+        zPostCenter = nbrZPost[0]
     }
 
-    if (lrRaw <= 0 || lrRaw >= maxRatio) return hw.lrRatioInvalid
+    // 6 boundary points: 2 center (x=0) + 4 lateral
+    val xw = (listOf(0.0, 0.0) + latPtsX).toDoubleArray()
+    val zw = (listOf(zAntCenter, zPostCenter) + latPtsZ).toDoubleArray()
+
+    // Normalized ellipse fit: α·x̂² + β·ẑ² + γ·x̂ + δ·ẑ = 1
+    val xm = xw.average(); val xs = std(xw) + 1e-12
+    val zm = zw.average(); val zs = std(zw) + 1e-12
+    val xn = DoubleArray(xw.size) { (xw[it] - xm) / xs }
+    val zn = DoubleArray(zw.size) { (zw[it] - zm) / zs }
+
+    val sol = solveEllipseLSQ(xn, zn, xw.size) ?: return hw.lrRatioInvalid
+    val (alphaF, betaF, _, _) = sol.let { Triple(it[0], it[1], Pair(it[2], it[3])) }
+        .let { doubleArrayOf(sol[0], sol[1], sol[2], sol[3]) }
+
+    if (sol[0] <= 1e-12 || sol[1] <= 1e-12) return hw.lrRatioInvalid
+    val rVal = 1.0 + sol[2] * sol[2] / (4.0 * sol[0]) + sol[3] * sol[3] / (4.0 * sol[1])
+    if (rVal <= 0) return hw.lrRatioInvalid
+
+    val bLr = sqrt(rVal / sol[0]) * xs   // LR 반축 (mm)
+    val aAp = sqrt(rVal / sol[1]) * zs   // AP 반축 (mm)
+    if (aAp <= 0) return hw.lrRatioInvalid
+
+    val lrRaw = max(bLr / aAp, 1.0)
 
     // Shrinkage toward prior
-    val avgRatio = latResults.map { it.ratio }.average()
+    val dCenterMean = abs(zPostCenter - zAntCenter)
+    if (dCenterMean <= 0) return hw.lrRatioInvalid
+    val avgRatio = if (latChords.isNotEmpty()) latChords.average() / dCenterMean else 1.0
     val confidence = ((1.0 - avgRatio) / 0.10).coerceIn(0.0, 1.0)
     val result = hw.lrPrior + (lrRaw - hw.lrPrior) * confidence
-
-    // 하한 1.0: 원형보다 좁은 단면은 물리적으로 불가
     return result.coerceAtLeast(1.0)
 }
 
@@ -266,70 +315,6 @@ private fun segmentToDistancesMm(
     return Pair(min(d1, d2), max(d1, d2))
 }
 
-// ── Ellipse Cap Height Fitting ──
-// bv_estimation.py _ellipse_cap_heights() 1:1 포팅
-// 8개 경계점(4ch × ant/post)으로 축 정렬 타원 피팅 → cap 높이
-
-private data class EllipseCapResult(
-    val hBot: Double,
-    val hTop: Double,
-    val botKind: String,
-    val topKind: String
-)
-
-private fun ellipseCapHeights(
-    dAnt: List<Double>, dPost: List<Double>,
-    validChannels: List<Int>,
-    sensorZMm: DoubleArray, degreeDeg: DoubleArray,
-    aCapBot: Double, aCapTop: Double,
-    yS: DoubleArray
-): EllipseCapResult {
-    val ptsX = mutableListOf<Double>()
-    val ptsY = mutableListOf<Double>()
-
-    for ((i, ch) in validChannels.withIndex()) {
-        val theta = degreeDeg[ch] * Math.PI / 180.0
-        ptsX.add(dAnt[i] * cos(theta))
-        ptsY.add(sensorZMm[ch] + dAnt[i] * sin(theta))
-        ptsX.add(dPost[i] * cos(theta))
-        ptsY.add(sensorZMm[ch] + dPost[i] * sin(theta))
-    }
-
-    if (ptsX.size < 5) {
-        return EllipseCapResult(aCapBot, aCapTop, "fallback", "fallback")
-    }
-
-    // 축 정렬 타원: Ax² + Cy² + Dx + Ey = 1 (F=-1 정규화)
-    // least squares: M @ [A,C,D,E]^T = 1
-    val nPts = ptsX.size
-    val sol = solveEllipseLSQ(ptsX.toDoubleArray(), ptsY.toDoubleArray(), nPts)
-        ?: return EllipseCapResult(aCapBot, aCapTop, "fallback", "fallback")
-
-    val (aa, cc, dd, ee) = sol
-    if (aa <= 0 || cc <= 0) {
-        return EllipseCapResult(aCapBot, aCapTop, "fallback", "fallback")
-    }
-
-    val xc = -dd / (2 * aa)
-    val yc = -ee / (2 * cc)
-    val rhs = dd * dd / (4 * aa) + ee * ee / (4 * cc) + 1.0
-    if (rhs <= 0) {
-        return EllipseCapResult(aCapBot, aCapTop, "fallback", "fallback")
-    }
-
-    val cSemi = sqrt(rhs / cc) // SI 반축
-    val bladderBot = yc - cSemi
-    val bladderTop = yc + cSemi
-
-    var hBot = yS[0] - bladderBot
-    var hTop = bladderTop - yS[yS.size - 1]
-
-    hBot = hBot.coerceIn(0.0, aCapBot)
-    hTop = hTop.coerceIn(0.0, aCapTop)
-
-    return EllipseCapResult(hBot, hTop, "ellipse", "ellipse")
-}
-
 /** 4×4 least squares: M^T M x = M^T 1 */
 private fun solveEllipseLSQ(px: DoubleArray, py: DoubleArray, n: Int): DoubleArray? {
     // Build 4×4 normal equations: (M^T M) params = M^T ones
@@ -390,9 +375,10 @@ private fun solve4x4(A: Array<DoubleArray>, b: DoubleArray): DoubleArray? {
 fun estimateBladderVolume6ch(
     allWalls: List<Pair<Int, Int>?>,   // 6채널 전체 (ant, post), null = invalid
 ): BVResult? {
-    val centerWalls = PiezoHW.centerCh.map { if (it < allWalls.size) allWalls[it] else null }
+    var centerWalls = PiezoHW.centerCh.map { if (it < allWalls.size) allWalls[it] else null }
+    centerWalls = repairCenterWallsFor6ch(centerWalls)
     val lateralWalls = PiezoHW.lateralCh.map { if (it < allWalls.size) allWalls[it] else null }
-    val lrRatio = computeLrRatio(centerWalls, lateralWalls)
+    val lrRatio = computeLrRatio(centerWalls, lateralWalls, sensorZMm = PiezoHW.sensorZMmAll)
     return estimateBladderVolume(
         walls = centerWalls,
         lrRatio = lrRatio,
@@ -472,7 +458,24 @@ fun estimateBladderVolume(
         }
         return null
     }
-    val n = validChannels.size
+    var n = validChannels.size
+
+    // 2-1) Post median outlier 제거 (#21 merge)
+    if (n >= 3) {
+        val posts = dPost.toDoubleArray()
+        val med = posts.sorted()[posts.size / 2]
+        val postTol = max(med * 0.25, 5.0 * distancePerSample)
+        val keep = (0 until n).filter { abs(posts[it] - med) <= postTol }
+        if (keep.size >= 2 && keep.size < n) {
+            val newValid = keep.map { validChannels[it] }.toMutableList()
+            val newDAnt = keep.map { dAnt[it] }.toMutableList()
+            val newDPost = keep.map { dPost[it] }.toMutableList()
+            validChannels.clear(); validChannels.addAll(newValid)
+            dAnt.clear(); dAnt.addAll(newDAnt)
+            dPost.clear(); dPost.addAll(newDPost)
+            n = validChannels.size
+        }
+    }
 
     // theta, sensor_z for valid channels
     val theta = DoubleArray(n) { degreeDeg[validChannels[it]] * Math.PI / 180.0 }
@@ -486,23 +489,70 @@ fun estimateBladderVolume(
     val S = DoubleArray(n) { areaK * D[it] * D[it] * lrRatio }
     val aCap = DoubleArray(n) { D[it] / 2.0 }
 
-    // 5) Midpoint 좌표
+    // 5) Midpoint + 벽 좌표 (y-z 평면, 타원 피팅용)
     val dMid = DoubleArray(n) { (dAnt[it] + dPost[it]) / 2.0 }
     val y = DoubleArray(n) { sensZ[it] + dMid[it] * sin(theta[it]) }
+    val yWallAnt = DoubleArray(n) { sensZ[it] + dAnt[it] * sin(theta[it]) }
+    val zWallAnt = DoubleArray(n) { dAnt[it] * cos(theta[it]) }
+    val yWallPost = DoubleArray(n) { sensZ[it] + dPost[it] * sin(theta[it]) }
+    val zWallPost = DoubleArray(n) { dPost[it] * cos(theta[it]) }
 
     // 6) y 오름차순 정렬
     val order = (0 until n).sortedBy { y[it] }
     val yS = DoubleArray(order.size) { y[order[it]] }
-    var sS = DoubleArray(order.size) { S[order[it]] }
-    var aCapS = DoubleArray(order.size) { aCap[order[it]] }
+    val sS = DoubleArray(order.size) { S[order[it]] }
+    val aCapS = DoubleArray(order.size) { aCap[order[it]] }
     val sortedCh = order.map { validChannels[it] }
 
-    // 7) 타원 피팅 → cap 높이
-    val ellCap = ellipseCapHeights(
-        dAnt, dPost, validChannels,
-        sensorZMm, degreeDeg,
-        aCapS[0], aCapS[n - 1], yS
-    )
+    // 7) y-z 평면 타원 피팅 → cap 높이 추정 (#21 merge)
+    val yw = DoubleArray(2 * n) { i -> if (i < n) yWallAnt[i] else yWallPost[i - n] }
+    val zw = DoubleArray(2 * n) { i -> if (i < n) zWallAnt[i] else zWallPost[i - n] }
+    val nPts = yw.size
+    var capKind = "fallback"
+    var z0Ellipse: Double? = null
+    var hCapBot = aCapS[0]
+    var hCapTop = aCapS[n - 1]
+
+    if (nPts >= 4) {
+        val ym = yw.average(); val ysStd = std(yw) + 1e-12
+        val zm = zw.average(); val zsStd = std(zw) + 1e-12
+        val yn = DoubleArray(nPts) { (yw[it] - ym) / ysStd }
+        val zn = DoubleArray(nPts) { (zw[it] - zm) / zsStd }
+
+        val sol = solveEllipseLSQ(yn, zn, nPts)
+        if (sol != null && sol[0] > 1e-12 && sol[1] > 1e-12) {
+            val y0n = -sol[2] / (2.0 * sol[0])
+            val z0n = -sol[3] / (2.0 * sol[1])
+            val rVal = 1.0 + sol[2] * sol[2] / (4.0 * sol[0]) + sol[3] * sol[3] / (4.0 * sol[1])
+            if (rVal > 0) {
+                val bSiN = sqrt(rVal / sol[0])
+                val y0 = y0n * ysStd + ym
+                z0Ellipse = z0n * zsStd + zm
+                val bSi = bSiN * ysStd
+                hCapBot = max(0.0, yS[0] - (y0 - bSi))
+                hCapTop = max(0.0, (y0 + bSi) - yS[n - 1])
+                hCapBot = min(hCapBot, aCapS[0])
+                hCapTop = min(hCapTop, aCapS[n - 1])
+                capKind = "ellipse"
+            }
+        }
+    }
+
+    // Top cap 상한: 미검출 상위 채널의 빔 y 좌표로 제한
+    val nTotalCh = sensorZMm.size
+    val topSortedCh = sortedCh.last()
+    if (n < nTotalCh && topSortedCh > 0) {
+        val upperCh = topSortedCh - 1
+        val upperTheta = degreeDeg[upperCh] * Math.PI / 180.0
+        val yBeamUpper = if (z0Ellipse != null && abs(cos(upperTheta)) > 1e-6) {
+            val dAtZ0 = z0Ellipse!! / cos(upperTheta)
+            sensorZMm[upperCh] + dAtZ0 * sin(upperTheta)
+        } else {
+            sensorZMm[upperCh].toDouble()
+        }
+        val hTopLimit = yBeamUpper - yS[n - 1]
+        if (hTopLimit > 0 && hCapTop > hTopLimit) hCapTop = hTopLimit
+    }
 
     // 8) Core frustum (traditional: h = dy)
     val dy = DoubleArray(n - 1) { yS[it + 1] - yS[it] }
@@ -511,9 +561,11 @@ fun estimateBladderVolume(
     }
     val vCore = vFrustum.sum()
 
-    // 9) Caps — Bottom: sphere, Top: cone
-    val vBottom = sS[0] * ellCap.hBot / 2.0 + Math.PI * ellCap.hBot * ellCap.hBot * ellCap.hBot / 6.0
-    val vTop = sS[n - 1] * ellCap.hTop / 3.0
+    // 9) Caps — Bottom: spherical cap, Top: cone
+    val vBottom = sS[0] * hCapBot / 2.0 + Math.PI * hCapBot * hCapBot * hCapBot / 6.0
+    val vTop = sS[n - 1] * hCapTop / 3.0
+    val bottomKind = "$capKind sphere"
+    val topKind = "$capKind cone"
 
     // 10) 합산
     val bvMm3 = vCore + vBottom + vTop
@@ -533,10 +585,10 @@ fun estimateBladderVolume(
         vCoreMm3 = vCore,
         vBottomMm3 = vBottom,
         vTopMm3 = vTop,
-        bottomHMm = ellCap.hBot,
-        topHMm = ellCap.hTop,
-        bottomKind = ellCap.botKind,
-        topKind = ellCap.topKind,
+        bottomHMm = hCapBot,
+        topHMm = hCapTop,
+        bottomKind = bottomKind,
+        topKind = topKind,
         lrRatio = lrRatio,
         distancePerSample = distancePerSample,
         delayOffsetMm = delayOffsetMm