From 139ebcbbb1852318fe8e5a6011108f0552eaa143 Mon Sep 17 00:00:00 2001 From: jjangddu Date: Thu, 7 May 2026 12:06:41 +0900 Subject: [PATCH] =?UTF-8?q?refactor:=20BV=20estimation=20=E2=80=94=20y-z/x?= =?UTF-8?q?-z=20ellipse=20fitting,=20wall=20repair,=20outlier=20removal?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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) --- .../managers/PiezoBVEstimator.kt | 340 ++++++++++-------- 1 file changed, 196 insertions(+), 144 deletions(-) diff --git a/app/src/main/java/com/example/medilightv2android/managers/PiezoBVEstimator.kt b/app/src/main/java/com/example/medilightv2android/managers/PiezoBVEstimator.kt index 2250146..0a5d643 100644 --- a/app/src/main/java/com/example/medilightv2android/managers/PiezoBVEstimator.kt +++ b/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?> +): List?> { + 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?>, lateralWalls: List?>, - 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() - 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() + // Lateral wall points (x-z plane) + val latPtsX = mutableListOf() + val latPtsZ = mutableListOf() + val latSiPositions = mutableListOf() + val latChords = mutableListOf() 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() + val nbrZAnt = mutableListOf() + val nbrZPost = mutableListOf() + 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, dPost: List, - validChannels: List, - sensorZMm: DoubleArray, degreeDeg: DoubleArray, - aCapBot: Double, aCapTop: Double, - yS: DoubleArray -): EllipseCapResult { - val ptsX = mutableListOf() - val ptsY = mutableListOf() - - 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, b: DoubleArray): DoubleArray? { fun estimateBladderVolume6ch( allWalls: List?>, // 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