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 556f8ab..e9aac09 100644 --- a/app/src/main/java/com/medithings/vesiscan/managers/GreenZoneConstants.kt +++ b/app/src/main/java/com/medithings/vesiscan/managers/GreenZoneConstants.kt @@ -30,11 +30,12 @@ enum class PlacementGuideMode { SIMPLE, BOUNDARY, SWEEP } * - METHOD_D 에 phantom-style lr floor (1.0) 를 추가하지 말 것. 사용자 결정 * 2026-06-30: "phantom QC 는 V41 로, 인체는 Python 1:1". */ -// 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 } +// METHOD_D_PHANTOM: 2026-08-11 신설. 순수 sphere 공식 (V = 4/3·π·r³). +// 방광을 완전한 구로 가정 · walls (ant/post) 로부터 지름 D 를 평균해 반지름 산출. +// 목적: demo-final 브랜치 = phantom 시연 전용 (실측 300ml phantom · 지름 ~83mm → +// ~299ml 로 정확 재현). 인체 (타원 방광) 에는 부적합. +// cloud-mvp = 인체 실사용 · METHOD_D (Halir + adaptive) 유지. +enum class BvMethod { FRUSTUM, V41, METHOD_D, METHOD_D_PHANTOM } /** Sensor alignment 알고리즘 — V1=기존 (computePlacementGuide), V2=신규 (alignment.py 포팅). */ enum class AlignmentAlgo { V1, V2 } @@ -52,11 +53,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). - // 2026-08-11 (rev): default → METHOD_D 로 복귀. - // METHOD_D_P (legacy Frustum+cap) 는 Halir-Flusser · adaptive_large_bladder_relax - // 미탑재라 큰 방광 (~300ml phantom) 에서 원리상 값이 150 대까지 낮게 나옴. - // 사용자 관찰 250ml 는 METHOD_D (Halir+adaptive 적용) 결과. dev 토글에는 D-P 남김. - @Volatile var bvMethod: BvMethod = BvMethod.METHOD_D + // 2026-08-11 (final): demo-final default = METHOD_D_PHANTOM (구 공식). + // Phantom 실측 300ml (지름 ~83mm) 시연에서 300ml 근처로 정확 재현. + // cloud-mvp 는 METHOD_D 유지 (인체 · Halir + adaptive) — 브랜치별로 다름이 정상 설계. + @Volatile var bvMethod: BvMethod = BvMethod.METHOD_D_PHANTOM /** * lr_ratio 강제 override (algorithm 팀 최신 표준: piezophantomtest 4830e7c). diff --git a/app/src/main/java/com/medithings/vesiscan/managers/PiezoBVEstimator.kt b/app/src/main/java/com/medithings/vesiscan/managers/PiezoBVEstimator.kt index 35e4f17..1faa9c2 100644 --- a/app/src/main/java/com/medithings/vesiscan/managers/PiezoBVEstimator.kt +++ b/app/src/main/java/com/medithings/vesiscan/managers/PiezoBVEstimator.kt @@ -583,6 +583,62 @@ fun estimateBladderVolume6ch( /** Wall + span 정보 (Python `extract_walls` 반환 4-tuple 대응). */ data class WallWithSpan(val ant: Double, val post: Double, val lowStart: Int, val lowEnd: Int) +/** + * 2026-08-11: demo-final (phantom 시연 전용) BV 계산. + * + * 방광을 **완전한 구 (sphere)** 로 가정. walls 의 (post-ant) 를 지름으로 보고 · 채널 + * 평균 지름 D 로부터 V = 4/3 · π · r³ 계산 (r = D/2). + * + * 목적: + * - phantom (실측 300 ml · 지름 ~83 mm) 시연에서 300 ml 근처 정확 재현. + * - cloud-mvp 인체용 (METHOD_D · Halir + adaptive) 와 격리된 데모 채널. + * + * 검증: + * - D = 42 samples × 1.968 dps = 82.7 mm → r = 41.35 mm → V ≈ 296 ml ✓ + * - 인체 (타원 방광) 에는 부적합. cloud-mvp 는 반드시 METHOD_D 사용. + */ +fun estimateBvPhantomSphere(walls: List): BVResult? { + val dps = PiezoHW.distancePerSample + // center 4 channels (CH0~CH3) 만 사용 · lateral (CH4/5) 제외. + val centerIdx: List = PiezoHW.centerCh.toList() + val validChannels: List = centerIdx.filter { idx -> + idx < walls.size && walls[idx] != null + } + val diameters: List = validChannels.map { idx -> + val w = walls[idx]!! + (w.post - w.ant) * dps + } + if (diameters.isEmpty()) return null + val d = diameters.average() + if (d <= 0.0) return null + val r = d / 2.0 + val volMm3 = (4.0 / 3.0) * Math.PI * r * r * r + val volMl = volMm3 / 1000.0 + return BVResult( + volumeMl = volMl, + volumeMm3 = volMm3, + validChannels = validChannels, + dAntMm = DoubleArray(0), + dPostMm = DoubleArray(0), + dMm = diameters.toDoubleArray(), + sMm2 = DoubleArray(0), + sortedChannels = validChannels, + sortedYMm = DoubleArray(0), + sortedSMm2 = DoubleArray(0), + vFrustumMm3 = DoubleArray(0), + vCoreMm3 = 0.0, + vBottomMm3 = 0.0, + vTopMm3 = 0.0, + bottomHMm = 0.0, + topHMm = 0.0, + bottomKind = "sphere", + topKind = "sphere", + lrRatio = 1.0, + distancePerSample = dps, + delayOffsetMm = PiezoHW.delayOffsetMm, + ) +} + /** * Python `runners.estimate_bv(walls, hw, adaptive_large_bladder_relax=True)` 1:1 이식. * 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 deleted file mode 100644 index 09ee72c..0000000 --- a/app/src/main/java/com/medithings/vesiscan/managers/legacydp/PiezoBVEstimatorLegacyDP.kt +++ /dev/null @@ -1,822 +0,0 @@ -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+ 기준). - // 2026-08-11: demo-final default 를 1.968 로 복귀 (phantom 시연 스케일 재현). - @Volatile private var _distancePerSample: Double = 1.968.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 6c9b6b4..5c2d913 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 @@ -427,9 +427,10 @@ fun PiezoMonitoringView(appState: AppState) { // 2026-06-30: BvMethod.METHOD_D 시 Method D walls 사용. 단 detection 이 // 이미 METHOD_D 면 allWalls 가 곧 Method D → recompute 생략. - // 2026-08-11: METHOD_D_P (legacy phantom fork) 도 동일 경로 통과 · dispatcher 에서 분기. + // 2026-08-11: METHOD_D_PHANTOM (구 공식) 도 동일 경로 통과 · dispatcher 에서 분기. + // Method D walls 를 그대로 sphere 계산에 사용 (D=post-ant 평균). val useMethodDBv = bvMethodSetting == com.medithings.vesiscan.managers.BvMethod.METHOD_D || - bvMethodSetting == com.medithings.vesiscan.managers.BvMethod.METHOD_D_P + bvMethodSetting == com.medithings.vesiscan.managers.BvMethod.METHOD_D_PHANTOM val detectionIsMethodD = method == com.medithings.vesiscan.managers.DetectionMethod.METHOD_D val sourceAllWalls: List?> = if (useMethodDBv && !detectionIsMethodD) { val signals = (0..5).map { ch -> @@ -483,15 +484,11 @@ fun PiezoMonitoringView(appState: AppState) { lowStart = it.lowStart, lowEnd = it.lowEnd, ) } } - // 2026-08-11: BvMethod dispatcher — METHOD_D_P (phantom · 구 가정) - // 는 bcbc7b4 시점 legacy 로직 (Halir-Flusser 등 Python parity fix 이전) 사용. - // Legacy 함수 시그니처는 List?> 이라 WallWithSpan → Pair 변환. + // 2026-08-11: BvMethod dispatcher — METHOD_D_PHANTOM (구 공식 · demo 시연) + // vs METHOD_D (Halir + adaptive · 인체용 · cloud-mvp default). 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) + com.medithings.vesiscan.managers.BvMethod.METHOD_D_PHANTOM) { + com.medithings.vesiscan.managers.estimateBvPhantomSphere(walls) } else { com.medithings.vesiscan.managers.estimateBv(walls) } @@ -2048,7 +2045,7 @@ fun PiezoMonitoringView(appState: AppState) { com.medithings.vesiscan.managers.BvMethod.FRUSTUM to "Frustum", com.medithings.vesiscan.managers.BvMethod.V41 to "V41", com.medithings.vesiscan.managers.BvMethod.METHOD_D to "D", - com.medithings.vesiscan.managers.BvMethod.METHOD_D_P to "D-P" + com.medithings.vesiscan.managers.BvMethod.METHOD_D_PHANTOM to "D-Ph" ) val bvColors = listOf(MlPrimary, Color(0xFFFF5722), Color(0xFF7C3AED), Color(0xFF00897B)) Row(horizontalArrangement = Arrangement.spacedBy(4.dp)) {