From fd7b5dad53c10ad99c1009d50b3f5a16832c7e9b Mon Sep 17 00:00:00 2001 From: jjangddu Date: Tue, 11 Aug 2026 14:41:48 +0900 Subject: [PATCH] =?UTF-8?q?feat(bv):=20METHOD=5FD=5FP=20legacy=20fork=20?= =?UTF-8?q?=C2=B7=20phantom=20=EC=8B=9C=EC=97=B0=EC=9A=A9=20(bcbc7b4=20?= =?UTF-8?q?=EC=8B=9C=EC=A0=90=20=EB=A1=9C=EC=A7=81=20=EA=B2=A9=EB=A6=AC)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 배경: · 8/4 default 스위치 (bvMethod V41 → METHOD_D) 순간 · 7/20 도입돼 있던 Python parity fix 4종 (Halir-Flusser ellipse fit · adaptive_large_bladder_relax · b_si_floor · subsample refined wall) 이 처음 활성화됨. · 사용자 관찰 "이전 시연 대비 volume 값이 조금 작음" 은 이 fix 들의 정확도 개선 (실제값에 가까워짐) 효과. 인체엔 좋으나 phantom (구 형태 방광 모형) 시연 재현성이 손상됨. 해결: · BvMethod 에 METHOD_D_P 추가 (P = phantom · 구 가정) · demo-final default 로. · bcbc7b4 (2026-07-02) 시점 PiezoBVEstimator 를 legacydp 서브패키지에 격리. · BVResult · WallWithSpan · GreenZoneConstants 등 최상위 shared 심볼은 최신 참조 (BVResult 신규 optional 필드는 default 로 자동 채워짐). · PiezoHW 는 legacy 파일 내부에 격리 (bcbc7b4 시점 V0/V1/V2 preset · 최신 R100/R200/R300 preset 변경과 무관하게 phantom 재현성 유지). 파일: · [신규] managers/legacydp/PiezoBVEstimatorLegacyDP.kt (821 lines) - bcbc7b4 원본 855 lines 에서 BVResult 정의 (40 lines) 제거 · import 로 대체 - package = com.medithings.vesiscan.managers.legacydp · [수정] managers/GreenZoneConstants.kt — BvMethod 확장 · demo default 변경 · [수정] ui/views/monitoring/PiezoMonitoringView.kt — dispatcher 분기 (WallWithSpan → Pair 변환 · legacy estimateBladderVolume6ch 호출) 효과: · demo-final default = METHOD_D_P → phantom 값 재현 · 인체용/개발용은 dev-mode 토글로 METHOD_D 선택 가능 (기존 로직 무손상) · cloud-mvp 브랜치의 METHOD_D 개선사항은 계속 원본 코드로 흘러올 수 있음 (legacy 는 완전 격리 · touch 안 됨). 빌드: BUILD SUCCESSFUL 8s (첫 시도 통과). --- .../vesiscan/managers/GreenZoneConstants.kt | 11 +- .../legacydp/PiezoBVEstimatorLegacyDP.kt | 821 ++++++++++++++++++ .../views/monitoring/PiezoMonitoringView.kt | 13 +- 3 files changed, 842 insertions(+), 3 deletions(-) create mode 100644 app/src/main/java/com/medithings/vesiscan/managers/legacydp/PiezoBVEstimatorLegacyDP.kt diff --git a/app/src/main/java/com/medithings/vesiscan/managers/GreenZoneConstants.kt b/app/src/main/java/com/medithings/vesiscan/managers/GreenZoneConstants.kt index a69c146..4a1a2bb 100644 --- a/app/src/main/java/com/medithings/vesiscan/managers/GreenZoneConstants.kt +++ b/app/src/main/java/com/medithings/vesiscan/managers/GreenZoneConstants.kt @@ -30,7 +30,11 @@ enum class PlacementGuideMode { SIMPLE, BOUNDARY, SWEEP } * - METHOD_D 에 phantom-style lr floor (1.0) 를 추가하지 말 것. 사용자 결정 * 2026-06-30: "phantom QC 는 V41 로, 인체는 Python 1:1". */ -enum class BvMethod { FRUSTUM, V41, METHOD_D } +// METHOD_D_P: 2026-08-11 신설. Python parity fix (Halir-Flusser ellipse fit · +// adaptive_large_bladder_relax · b_si_floor · subsample refined 등) 도입 이전 +// (bcbc7b4 · 2026-07-02) 시점 BV 로직 격리 fork. phantom (구 가정) 시연 재현용. +// demo-final default. 인체 (타원 방광) 정확도 필요 시 METHOD_D 사용. +enum class BvMethod { FRUSTUM, V41, METHOD_D, METHOD_D_P } /** Sensor alignment 알고리즘 — V1=기존 (computePlacementGuide), V2=신규 (alignment.py 포팅). */ enum class AlignmentAlgo { V1, V2 } @@ -48,7 +52,10 @@ object GreenZoneConstants { // detectionMethod = METHOD_C (walls 검출 · phantom-검증 파이프라인) // bvMethod = METHOD_D (BV 계산 · Python parity + adaptive_large_bladder_relax + b_si_floor) // Method D BV 는 estimateBv(walls) wrapper 로 라우팅 (PiezoMonitoringView). - @Volatile var bvMethod: BvMethod = BvMethod.METHOD_D + // 2026-08-11: demo-final default = METHOD_D_P (phantom · 구 가정 · bcbc7b4 로직). + // 사용자 관찰 "이전 시연 대비 volume 값이 조금 작음" 은 7/20 Python parity fix 들이 + // 8/4 default 스위치 순간 활성화된 결과. 시연 재현성 위해 legacy fork 격리. + @Volatile var bvMethod: BvMethod = BvMethod.METHOD_D_P /** * lr_ratio 강제 override (algorithm 팀 최신 표준: piezophantomtest 4830e7c). diff --git a/app/src/main/java/com/medithings/vesiscan/managers/legacydp/PiezoBVEstimatorLegacyDP.kt b/app/src/main/java/com/medithings/vesiscan/managers/legacydp/PiezoBVEstimatorLegacyDP.kt new file mode 100644 index 0000000..a515734 --- /dev/null +++ b/app/src/main/java/com/medithings/vesiscan/managers/legacydp/PiezoBVEstimatorLegacyDP.kt @@ -0,0 +1,821 @@ +package com.medithings.vesiscan.managers.legacydp + +import com.medithings.vesiscan.walldetect.core.WdConfig +// 2026-08-11: legacy fork for METHOD_D_P (Python parity fix 미적용 · phantom 구 가정) +// Base: demo-final bcbc7b4 (2026-07-02) — Halir-Flusser ellipse fit 도입 전 상태. +// BVResult · WallWithSpan 은 최신 참조 (신규 optional 필드 자동 default). +// PiezoHW 는 이 파일 내부에 격리 (bcbc7b4 시점 · phantom 시연 안정성). +import com.medithings.vesiscan.managers.BVResult +import com.medithings.vesiscan.managers.WallWithSpan +import kotlin.math.abs +import kotlin.math.cos +import kotlin.math.hypot +import kotlin.math.min +import kotlin.math.max +import kotlin.math.sin +import kotlin.math.sqrt + +// ========================================================================== +// Bladder Volume Estimation — Multi-Channel Frustum + Cap (5ch vertical) +// +// 1:1 port of PiezoBVEstimator.swift / bv_estimation.py (algo-test branch) +// Only "traditional" mode is used (coord/hybrid disabled in Python too) +// ========================================================================== + +// ── Hardware Parameters (5ch vertical array) ── + +/** + * Probe geometry — probe model별 고정 + */ +object PiezoHW { + // ═══════════════════════════════════════════════════════════ + // Device Presets — 6ch 기기 3가지 버전 (Snell's law 굴절 보정 적용) + // config_6ch.py 1:1 포팅 + // ═══════════════════════════════════════════════════════════ + + enum class DevicePreset { + LEGACY_5CH, + V0, // All 0° (BLE 테스트용) + V1, // Max 20° housing → Snell's law refraction + V2 // Max 30° housing → Snell's law refraction + } + + var activePreset: DevicePreset = DevicePreset.V0 + + fun autoDetectPreset(deviceName: String) { + if (deviceName.startsWith("VBT") && deviceName.length >= 4) { + // VBT...n0x → 뒤에서 3번째 글자(n)가 각도 타입 + // n=0 → V0, n=2 → V1, n=3 → V2 + val n = deviceName[deviceName.length - 3] + activePreset = when (n) { + '0' -> DevicePreset.V0 + '2' -> DevicePreset.V1 + '3' -> DevicePreset.V2 + else -> DevicePreset.V1 // 알 수 없으면 V1 기본 + } + } else if (deviceName.startsWith("2025MEDIP")) { + activePreset = DevicePreset.LEGACY_5CH + } + android.util.Log.d("PiezoHW", "autoDetectPreset: '$deviceName' → ${activePreset.name}") + } + + val centerCh = intArrayOf(0, 1, 2, 3) + val lateralCh = intArrayOf(4, 5) + val lateralNeighbors = intArrayOf(1, 2) + + // 6채널 전체 z 좌표 (CH0~CH5) + val sensorZMmAll: DoubleArray + get() = when (activePreset) { + DevicePreset.LEGACY_5CH -> doubleArrayOf(22.0, 16.5, 11.0, 5.5, 0.0, 0.0) + DevicePreset.V0 -> doubleArrayOf(21.0, 14.0, 7.0, 0.0, 10.5, 10.5) + DevicePreset.V1 -> doubleArrayOf(20.1, 12.8, 7.0, 0.0, 9.9, 9.9) + DevicePreset.V2 -> doubleArrayOf(19.3, 13.0, 6.7, 0.0, 9.85, 9.85) + } + + // 6채널 전체 x 좌표 + val sensorXMmAll: DoubleArray + get() = when (activePreset) { + DevicePreset.LEGACY_5CH -> doubleArrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0) + else -> doubleArrayOf(0.0, 0.0, 0.0, 0.0, -10.0, 10.0) + } + + // 6채널 전체 SI 빔 각도 (Snell's law 굴절 후) + val degreeAll: DoubleArray + get() = when (activePreset) { + DevicePreset.LEGACY_5CH -> doubleArrayOf(0.0, -2.2, -4.4, -6.6, -8.8, 0.0) + DevicePreset.V0 -> doubleArrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0) + DevicePreset.V1 -> doubleArrayOf(6.89, 0.0, -6.89, -13.66, -6.87, -6.87) + DevicePreset.V2 -> doubleArrayOf(0.0, -6.89, -13.66, -20.20, -6.87, -6.87) + } + + // 6채널 전체 LR 빔 각도 (Snell's law 굴절 후) + val degreeLRAll: DoubleArray + get() = when (activePreset) { + DevicePreset.LEGACY_5CH -> doubleArrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0) + DevicePreset.V0 -> doubleArrayOf(0.0, 0.0, 0.0, 0.0, 0.0, 0.0) + DevicePreset.V1 -> doubleArrayOf(0.0, 0.0, 0.0, 0.0, -3.42, 3.42) + DevicePreset.V2 -> doubleArrayOf(0.0, 0.0, 0.0, 0.0, -3.42, 3.42) + } + + // BV용 center 채널만 (CH0~CH3) + val sensorZMm: DoubleArray get() = centerCh.map { sensorZMmAll[it] }.toDoubleArray() + val sensorXMm: DoubleArray get() = centerCh.map { sensorXMmAll[it] }.toDoubleArray() + val degree: DoubleArray get() = centerCh.map { degreeAll[it] }.toDoubleArray() + val degreeLR: DoubleArray get() = centerCh.map { degreeLRAll[it] }.toDoubleArray() + + /** Acoustic calibration (런타임 조절 가능). + * Setter는 walldetect.core.WdConfig.DPS_DEFAULT 도 함께 동기화한다. + * V41Detector / Geometry / AnatomicalGate 가 WdConfig 측 상수를 사용하기 + * 때문에 한 곳만 바꾸면 BV 두 경로 (6ch estimate vs V41 dispatch / gate) 가 + * 서로 다른 dps 로 계산되어 측정이 어긋난다. */ + // `.also` 로 초기값을 walldetect.core.WdConfig 측에도 동시에 반영 — 시작 + // 시점부터 두 경로의 dps 가 일치한다. + // 02bed02 (config_6ch.py): DISTANCE_PER_SAMPLE 1.981 → 1.936 (HW 샘플링레이트 변경). + // 우리 default 도 1.968 → 1.936 으로 정렬 (FW VBTFW118+ 기준). + @Volatile private var _distancePerSample: Double = 1.936.also { WdConfig.DPS_DEFAULT = it } + var distancePerSample: Double + get() = _distancePerSample + set(value) { + _distancePerSample = value + WdConfig.DPS_DEFAULT = value + } + const val delayOffsetMm: Double = 6.85 + + /** Volume model */ + val areaK: Double = Math.PI / 4.0 + const val defaultLrRatio: Double = 1.0 + + // lr_ratio fallbacks + const val lrRatioNoDetection: Double = 1.0 + const val lrRatioInvalid: Double = 1.2 + const val lrPrior: Double = 1.2 + + val presetName: String + get() = when (activePreset) { + DevicePreset.LEGACY_5CH -> "Legacy 5ch" + DevicePreset.V0 -> "V0 (all 0°)" + DevicePreset.V1 -> "V1 (max 20°, Snell)" + DevicePreset.V2 -> "V2 (max 30°, Snell)" + } +} + +// ── 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, + sensorZMm: DoubleArray = PiezoHW.sensorZMmAll +): Double { + val hw = PiezoHW + val neighbors = hw.lateralNeighbors + val availableNbrs = neighbors.filter { ni -> + ni < centerWalls.size && centerWalls[ni] != null + } + if (availableNbrs.size < 2) return hw.lrRatioNoDetection + + // 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 + val w = lateralWalls[li] ?: continue + val (dNear, dFar) = segmentToDistancesMm(w.first, w.second) + 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() + + // 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)) + } + + 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 { + zAntCenter = nbrZAnt[0] + zPostCenter = nbrZPost[0] + } + + // 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 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 + return result.coerceAtLeast(1.0) +} + +// ── Result ── + + +// ── Helpers ── + +/** Sample index → distance (mm) */ +private fun sampleToMm( + idx: Double, + dps: Double = PiezoHW.distancePerSample, + offset: Double = PiezoHW.delayOffsetMm +): Double = offset + idx * dps + +/** (ant, post) sample indices → (d_near, d_far) mm */ +private fun segmentToDistancesMm( + ant: Int, post: Int, + dps: Double = PiezoHW.distancePerSample, + offset: Double = PiezoHW.delayOffsetMm +): Pair { + val d1 = sampleToMm(ant.toDouble(), dps, offset) + val d2 = sampleToMm(post.toDouble(), dps, offset) + return Pair(min(d1, d2), max(d1, d2)) +} + +/** 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 + // M columns: [x², y², x, y] + val mtm = Array(4) { DoubleArray(4) } + val mtb = DoubleArray(4) + + for (k in 0 until n) { + val x = px[k]; val y = py[k] + val row = doubleArrayOf(x * x, y * y, x, y) + for (i in 0 until 4) { + for (j in 0 until 4) mtm[i][j] += row[i] * row[j] + mtb[i] += row[i] // RHS = 1 + } + } + return solve4x4(mtm, mtb) +} + +private fun solve4x4(A: Array, b: DoubleArray): DoubleArray? { + val a = Array(4) { A[it].copyOf() } + val bb = b.copyOf() + for (col in 0 until 4) { + var maxRow = col; var maxVal = abs(a[col][col]) + for (row in (col + 1) until 4) { + if (abs(a[row][col]) > maxVal) { maxVal = abs(a[row][col]); maxRow = row } + } + if (maxVal < 1e-12) return null + if (maxRow != col) { + val tmpA = a[col]; a[col] = a[maxRow]; a[maxRow] = tmpA + val tmpB = bb[col]; bb[col] = bb[maxRow]; bb[maxRow] = tmpB + } + for (row in (col + 1) until 4) { + val factor = a[row][col] / a[col][col] + for (j in col until 4) a[row][j] -= factor * a[col][j] + bb[row] -= factor * bb[col] + } + } + val x = DoubleArray(4) + for (i in 3 downTo 0) { + var sum = bb[i] + for (j in (i + 1) until 4) sum -= a[i][j] * x[j] + if (abs(a[i][i]) < 1e-12) return null + x[i] = sum / a[i][i] + } + return x +} + +// ── NEW BV lib (piezophantomtest 4830e7c) helpers ───────────────────────── + +/** posterior arc chord-대비 최대 수직 이탈 (sagitta, mm). Python `_sagitta_mm`. */ +private fun sagittaMm(yp: DoubleArray, zp: DoubleArray): Double { + if (yp.size < 3) return 0.0 + val order = (0 until yp.size).sortedBy { yp[it] } + val y = DoubleArray(yp.size) { yp[order[it]] } + val z = DoubleArray(zp.size) { zp[order[it]] } + val p0y = y[0]; val p0z = z[0] + val dY = y[y.size - 1] - p0y + val dZ = z[z.size - 1] - p0z + val ln = hypot(dY, dZ) + if (ln < 1e-6) return 0.0 + val nY = -dZ / ln; val nZ = dY / ln + var m = 0.0 + for (i in 1 until y.size - 1) { + val d = abs((y[i] - p0y) * nY + (z[i] - p0z) * nZ) + if (d > m) m = d + } + return m +} + +/** + * posterior arc 원 fit (Kåsa algebraic) + sagitta shrinkage → R_eff. + * Python `_shrink_si_radius`. 실패 시 null. + */ +private fun shrinkSiRadius( + yWallPost: DoubleArray, zWallPost: DoubleArray, + cApPrior: Double?, + sNoise: Double = BOTTOM_CAP_SAGITTA_NOISE_MM, +): Double? { + val n = yWallPost.size + if (n < 3 || cApPrior == null || cApPrior <= 0) return null + val ata = Array(3) { DoubleArray(3) } + val atb = DoubleArray(3) + for (i in 0 until n) { + val row = doubleArrayOf(yWallPost[i], zWallPost[i], 1.0) + val r = yWallPost[i] * yWallPost[i] + zWallPost[i] * zWallPost[i] + for (a in 0 until 3) { + for (b in 0 until 3) ata[a][b] += row[a] * row[b] + atb[a] += row[a] * r + } + } + val sol = solve3x3(ata, atb) ?: return null + val yc = sol[0] / 2.0 + val zc = sol[1] / 2.0 + val rFit = sqrt(max(sol[2] + yc * yc + zc * zc, 1e-9)) + val s = sagittaMm(yWallPost, zWallPost) + val w = (s * s) / (s * s + sNoise * sNoise) + val kappa = w / rFit + (1.0 - w) / cApPrior + return if (kappa > 1e-9) 1.0 / kappa else cApPrior +} + +/** R 구에서 base a 인 minor 구면 캡 높이. Python `_minor_cap_height`. */ +private fun minorCapHeight(rEff: Double, aBase: Double): Double { + val a = min(aBase, rEff) + return rEff - sqrt(max(rEff * rEff - a * a, 0.0)) +} + +/** 구면 캡 부피 (lr 보정). V = π·h²·(3R−h)/3 · lr. Python `_spherical_cap_volume`. */ +private fun sphericalCapVolume(rEff: Double, h: Double, lrRatio: Double): Double = + Math.PI * h * h * (3.0 * rEff - h) / 3.0 * lrRatio + +/** + * 벽 인셋: (ant + f·(ls-ant), post - f·(post-le)) → 벽을 lumen 내부로 살짝 밀어 넣기. + * Python `_apply_lumen_inset_one`. frac=0 이면 (ant, post) 그대로. + * ls, le 없으면 (ant, post) 그대로. + */ +fun applyLumenInsetOne( + ant: Int, post: Int, ls: Int?, le: Int?, frac: Double = LUMEN_INSET_FRAC, +): Pair { + if (frac <= 0 || ls == null || le == null) return Pair(ant, post) + val ai = kotlin.math.round(ant + frac * (ls - ant)).toInt() + val pi = kotlin.math.round(post - frac * (post - le)).toInt() + return Pair(ai, pi) +} + +/** Python `BOTTOM_CAP_SAGITTA_NOISE_MM` (bv_estimation.py:895). */ +const val BOTTOM_CAP_SAGITTA_NOISE_MM: Double = 4.0 + +/** Python `LUMEN_INSET_FRAC` (config_6ch.py:41). */ +const val LUMEN_INSET_FRAC: Double = 0.15 + +// ── Core BV Computation ── + +/** + * Traditional frustum BV 추정 (cos 보정 + SI 축 기반) + * 1:1 port of _bv_core(mode="traditional") + */ +/** + * 6ch BV 추정 — center 4채널(CH0~CH3) + lateral 2채널(CH4/CH5)에서 lr_ratio 계산 + * config_6ch.py / bv_estimation.py estimate_bladder_volume_6ch 1:1 포팅 + */ +fun estimateBladderVolume6ch( + allWalls: List?>, // 6채널 전체 (ant, post), null = invalid +): BVResult? { + 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 } + // 2026-07-01: algorithm 팀 최신 표준 (piezophantomtest 4830e7c) — lr_ratio_override=1.0 + // Phantom 검증: computed lr 은 6채널 HW 한계로 -38% 오차. Dev panel 로 전환 가능. + val lrRatio = com.medithings.vesiscan.managers.GreenZoneConstants.lrRatioOverride + ?: computeLrRatio(centerWalls, lateralWalls, sensorZMm = PiezoHW.sensorZMmAll) + return estimateBladderVolume( + walls = centerWalls, + lrRatio = lrRatio, + edgeCh = 3, + edgeRefChannels = listOf(1, 2) + ) +} + +fun estimateBladderVolume( + walls: List?>, + distancePerSample: Double = PiezoHW.distancePerSample, + delayOffsetMm: Double = PiezoHW.delayOffsetMm, + sensorZMm: DoubleArray = PiezoHW.sensorZMm, + degreeDeg: DoubleArray = PiezoHW.degree, + areaK: Double = PiezoHW.areaK, + lrRatio: Double = PiezoHW.defaultLrRatio, + applyAngleCorrection: Boolean = true, + edgeCh: Int = walls.size - 1, + edgeRefChannels: List = listOf(walls.size - 2, walls.size - 3) +): BVResult? { + + // 1) Edge channel filter — cap 방식 결정용 플래그 + var edgeIsShort = false + if (edgeCh < walls.size) { + val wEdge = walls[edgeCh] + if (wEdge != null) { + val lEdge = wEdge.second - wEdge.first + val refLens = mutableListOf() + for (ci in edgeRefChannels) { + if (ci in walls.indices) { + val w = walls[ci] + if (w != null) refLens.add(w.second - w.first) + } + } + val minRef = refLens.minOrNull() + if (minRef != null && lEdge.toDouble() < 0.9 * minRef.toDouble()) { + edgeIsShort = true + } + } + } + + // 2) 유효 채널 추출 + sample → mm + val validChannels = mutableListOf() + val dAnt = mutableListOf() + val dPost = mutableListOf() + + for ((i, w) in walls.withIndex()) { + if (w == null) continue + validChannels.add(i) + val (dNear, dFar) = segmentToDistancesMm(w.first, w.second, distancePerSample, delayOffsetMm) + dAnt.add(dNear) + dPost.add(dFar) + } + + // 1채널 케이스: 구(sphere) 부피 (Python _single_channel_bv) + if (validChannels.size < 2) { + if (validChannels.size == 1) { + val ch = validChannels[0] + val theta = degreeDeg[ch] * Math.PI / 180.0 + val dRaw = dPost[0] - dAnt[0] + val D = if (applyAngleCorrection) dRaw * cos(theta) else dRaw + val R = D / 2.0 + val bvMm3 = (4.0 / 3.0) * Math.PI * R * R * R + val S = areaK * D * D * lrRatio + val yMid = sensorZMm[ch] + (dAnt[0] + dPost[0]) / 2.0 * sin(theta) + return BVResult( + volumeMl = bvMm3 / 1000.0, volumeMm3 = bvMm3, + validChannels = validChannels, dAntMm = dAnt.toDoubleArray(), dPostMm = dPost.toDoubleArray(), + dMm = doubleArrayOf(D), sMm2 = doubleArrayOf(S), + sortedChannels = listOf(ch), sortedYMm = doubleArrayOf(yMid), sortedSMm2 = doubleArrayOf(S), + vFrustumMm3 = doubleArrayOf(), vCoreMm3 = 0.0, + vBottomMm3 = bvMm3, vTopMm3 = 0.0, + bottomHMm = R, topHMm = 0.0, + bottomKind = "sphere", topKind = "sphere", + lrRatio = lrRatio, distancePerSample = distancePerSample, delayOffsetMm = delayOffsetMm + ) + } + return null + } + 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 } + val sensZ = DoubleArray(n) { sensorZMm[validChannels[it]] } + + // 3) 단면 직경 (traditional: cos 보정) + val lRaw = DoubleArray(n) { dPost[it] - dAnt[it] } + val D = DoubleArray(n) { if (applyAngleCorrection) lRaw[it] * cos(theta[it]) else lRaw[it] } + + // 4) 단면적 + AP 반지름 + val S = DoubleArray(n) { areaK * D[it] * D[it] * lrRatio } + val aCap = DoubleArray(n) { D[it] / 2.0 } + + // 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]] } + val sS = DoubleArray(order.size) { S[order[it]] } + val aCapS = DoubleArray(order.size) { aCap[order[it]] } + val sortedCh = order.map { validChannels[it] } + + // 7) y-z 평면 타원 피팅 + 반복 outlier 제거 (#22 merge) + val allYw = DoubleArray(2 * n) { i -> if (i < n) yWallAnt[i] else yWallPost[i - n] } + val allZw = DoubleArray(2 * n) { i -> if (i < n) zWallAnt[i] else zWallPost[i - n] } + val nPts = allYw.size + var capKind = "fallback" + var z0Ellipse: Double? = null + var cApPrior: Double? = null // shrink_bottom_cap 용 AP 반축 + var hCapBot = aCapS[0] + var hCapTop = aCapS[n - 1] + val ellipseCostThr = 0.5 // 점당 평균 잔차 임계 + + // 타원 피팅 helper: 성공 시 (y0, z0, bSi, aAp, residuals) 반환 + fun fitEllipsePts(ywF: DoubleArray, zwF: DoubleArray): Array? { + val cnt = ywF.size + val ym = ywF.average(); val ysS = std(ywF) + 1e-12 + val zm = zwF.average(); val zsS = std(zwF) + 1e-12 + val yn = DoubleArray(cnt) { (ywF[it] - ym) / ysS } + val zn = DoubleArray(cnt) { (zwF[it] - zm) / zsS } + val sol = solveEllipseLSQ(yn, zn, cnt) ?: return null + if (sol[0] <= 1e-12 || sol[1] <= 1e-12) return null + val rr = 1.0 + sol[2] * sol[2] / (4.0 * sol[0]) + sol[3] * sol[3] / (4.0 * sol[1]) + if (rr <= 0) return null + val y0 = (-sol[2] / (2.0 * sol[0])) * ysS + ym + val z0 = (-sol[3] / (2.0 * sol[1])) * zsS + zm + val bSi = sqrt(rr / sol[0]) * ysS + val aAp = sqrt(rr / sol[1]) * zsS + val resid = DoubleArray(cnt) { + abs(((ywF[it] - y0) / bSi) * ((ywF[it] - y0) / bSi) + + ((zwF[it] - z0) / aAp) * ((zwF[it] - z0) / aAp) - 1.0) + } + return arrayOf(y0, z0, bSi, aAp, resid) + } + + if (nPts >= 5) { + val keep = BooleanArray(nPts) { true } + var fit = fitEllipsePts(allYw, allZw) + + if (fit != null) { + var y0 = fit[0] as Double; var z0 = fit[1] as Double + var bSi = fit[2] as Double; var aAp = fit[3] as Double + var resid = fit[4] as DoubleArray + var meanRes = resid.average() + + // 반복 outlier 제거: worst 점 하나씩, 최소 5점 유지 + while (meanRes > ellipseCostThr && keep.count { it } > 5) { + val worstLocal = resid.indices.maxByOrNull { resid[it] } ?: break + val activeIndices = keep.indices.filter { keep[it] } + keep[activeIndices[worstLocal]] = false + val keptY = keep.indices.filter { keep[it] }.map { allYw[it] }.toDoubleArray() + val keptZ = keep.indices.filter { keep[it] }.map { allZw[it] }.toDoubleArray() + val fit2 = fitEllipsePts(keptY, keptZ) ?: break + val meanRes2 = (fit2[4] as DoubleArray).average() + if (meanRes2 < meanRes) { + y0 = fit2[0] as Double; z0 = fit2[1] as Double + bSi = fit2[2] as Double; aAp = fit2[3] as Double + resid = fit2[4] as DoubleArray; meanRes = meanRes2 + } else break + } + + // 품질 판정: 평균 잔차 ≤ threshold + if (meanRes <= ellipseCostThr) { + z0Ellipse = z0 + // b_si 상한: a_ap × 1.3 (해부학적 SI/AP 비율 제한) + if (bSi > aAp * 1.3) bSi = aAp * 1.3 + 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]) + cApPrior = aAp + capKind = "ellipse" + } + } + } + + // shrink_bottom_cap=true (Python default) — posterior arc 원 fit + sagitta shrink 로 + // R_eff 공유 곡률 산출. 성공 시 both caps 높이를 정규화 (자유 b_si 제거). + val rEffCap: Double? = + if (cApPrior != null) shrinkSiRadius(yWallPost, zWallPost, cApPrior) else null + if (rEffCap != null) { + hCapBot = minorCapHeight(rEffCap, aCapS[0]) + hCapTop = minorCapHeight(rEffCap, aCapS[n - 1]) + } + + // 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] } + val vFrustum = DoubleArray(n - 1) { + (dy[it] / 3.0) * (sS[it] + sS[it + 1] + sqrt(sS[it] * sS[it + 1])) + } + val vCore = vFrustum.sum() + + // 9) Caps — Bottom: spherical cap (shrink 시 R_eff 공유), Top: cone + val vBottom: Double + val vTop: Double + val bottomKind: String + val topKind: String + if (rEffCap != null) { + vBottom = sphericalCapVolume(rEffCap, hCapBot, lrRatio) + vTop = sS[n - 1] * hCapTop / 3.0 + bottomKind = "shrink sphere" + topKind = "shrink cone" + } else { + // fallback (ellipse fit 실패): 기존 hemisphere 공식. + vBottom = sS[0] * hCapBot / 2.0 + Math.PI * hCapBot * hCapBot * hCapBot / 6.0 + vTop = sS[n - 1] * hCapTop / 3.0 + bottomKind = "$capKind sphere" + topKind = "$capKind cone" + } + + // 10) 합산 + val bvMm3 = vCore + vBottom + vTop + + return BVResult( + volumeMl = bvMm3 / 1000.0, + volumeMm3 = bvMm3, + validChannels = validChannels, + dAntMm = dAnt.toDoubleArray(), + dPostMm = dPost.toDoubleArray(), + dMm = D, + sMm2 = S, + sortedChannels = sortedCh, + sortedYMm = yS, + sortedSMm2 = sS, + vFrustumMm3 = vFrustum, + vCoreMm3 = vCore, + vBottomMm3 = vBottom, + vTopMm3 = vTop, + bottomHMm = hCapBot, + topHMm = hCapTop, + bottomKind = bottomKind, + topKind = topKind, + lrRatio = lrRatio, + distancePerSample = distancePerSample, + delayOffsetMm = delayOffsetMm + ) +} + +// ── Math Utilities ── + +/** + * Degree 2 polynomial fit: returns Triple(c2, c1, c0) where f(x) = c2*x² + c1*x + c0 + * Least squares via normal equations (Vandermonde) + */ +private fun polyfit2(x: DoubleArray, y: DoubleArray): Triple { + val n = x.size + if (n < 3) return Triple(0.0, 0.0, y.firstOrNull() ?: 0.0) + + // Build normal equations for Ax = b where A is Vandermonde [x^0, x^1, x^2] + val sx = DoubleArray(5) + val sy = DoubleArray(3) + for (i in 0 until n) { + val xi = x[i] + val yi = y[i] + var xp = 1.0 + for (k in 0 until 5) { + sx[k] += xp + xp *= xi + } + sy[0] += yi + sy[1] += yi * xi + sy[2] += yi * xi * xi + } + + // 3x3 system + val a = arrayOf( + doubleArrayOf(sx[0], sx[1], sx[2]), + doubleArrayOf(sx[1], sx[2], sx[3]), + doubleArrayOf(sx[2], sx[3], sx[4]) + ) + val b = doubleArrayOf(sy[0], sy[1], sy[2]) + + val sol = solve3x3(a, b) ?: return Triple(0.0, 0.0, y.firstOrNull() ?: 0.0) + return Triple(sol[2], sol[1], sol[0]) // (c2, c1, c0) +} + +/** Evaluate degree 2 polynomial: c2*x² + c1*x + c0 */ +private fun polyval2(c: Triple, x: Double): Double = + c.first * x * x + c.second * x + c.third + +/** Solve 3x3 linear system via Gaussian elimination with partial pivoting */ +private fun solve3x3(A: Array, b: DoubleArray): DoubleArray? { + val a = Array(3) { A[it].copyOf() } + val bb = b.copyOf() + + for (col in 0 until 3) { + // Partial pivoting + var maxRow = col + var maxVal = abs(a[col][col]) + for (row in (col + 1) until 3) { + if (abs(a[row][col]) > maxVal) { + maxVal = abs(a[row][col]) + maxRow = row + } + } + if (maxVal < 1e-12) return null + if (maxRow != col) { + val tmpA = a[col]; a[col] = a[maxRow]; a[maxRow] = tmpA + val tmpB = bb[col]; bb[col] = bb[maxRow]; bb[maxRow] = tmpB + } + + // Eliminate + for (row in (col + 1) until 3) { + val factor = a[row][col] / a[col][col] + for (j in col until 3) a[row][j] -= factor * a[col][j] + bb[row] -= factor * bb[col] + } + } + + // Back substitution + val x = DoubleArray(3) + for (i in 2 downTo 0) { + var sum = bb[i] + for (j in (i + 1) until 3) sum -= a[i][j] * x[j] + if (abs(a[i][i]) < 1e-12) return null + x[i] = sum / a[i][i] + } + return x +} + +/** Standard deviation */ +private fun std(arr: DoubleArray): Double { + val n = arr.size.toDouble() + if (n <= 0) return 0.0 + val mean = arr.sum() / n + val variance = arr.sumOf { (it - mean).let { d -> d * d } } / n + return sqrt(variance) +} diff --git a/app/src/main/java/com/medithings/vesiscan/ui/views/monitoring/PiezoMonitoringView.kt b/app/src/main/java/com/medithings/vesiscan/ui/views/monitoring/PiezoMonitoringView.kt index 3d27181..c313183 100644 --- a/app/src/main/java/com/medithings/vesiscan/ui/views/monitoring/PiezoMonitoringView.kt +++ b/app/src/main/java/com/medithings/vesiscan/ui/views/monitoring/PiezoMonitoringView.kt @@ -481,7 +481,18 @@ fun PiezoMonitoringView(appState: AppState) { lowStart = it.lowStart, lowEnd = it.lowEnd, ) } } - val bvResult = com.medithings.vesiscan.managers.estimateBv(walls) + // 2026-08-11: BvMethod dispatcher — METHOD_D_P (phantom · 구 가정) + // 는 bcbc7b4 시점 legacy 로직 (Halir-Flusser 등 Python parity fix 이전) 사용. + // Legacy 함수 시그니처는 List?> 이라 WallWithSpan → Pair 변환. + val bvResult = if (com.medithings.vesiscan.managers.GreenZoneConstants.bvMethod == + com.medithings.vesiscan.managers.BvMethod.METHOD_D_P) { + val pairs: List?> = walls.map { w -> + w?.let { Pair(it.ant.toInt(), it.post.toInt()) } + } + com.medithings.vesiscan.managers.legacydp.estimateBladderVolume6ch(pairs) + } else { + com.medithings.vesiscan.managers.estimateBv(walls) + } if (bvResult != null) { rawVolumeMl = bvResult.volumeMl rawLrRatio = bvResult.lrRatio