feat: Otsu adaptive threshold, Method A 포팅, BV ellipse cap, UI 정리
Algorithm: - Otsu 1D adaptive threshold 포팅 (span_utils.py otsu_1d, 기본 ON) - Method A 벽 검출 포팅 (PiezoEchoAnalyzerA.kt: plateau score + prominence) - Method A cross-validation (채널간 중앙값/gradient 보정) - BV cap 높이: parabolic → ellipse 피팅 (8점 축정렬 타원) - Bottom cap: sphere, Top cap: cone (Python bv_estimation.py 동기화) - 설정에 Method A/B 토글 (왼쪽=A, 오른쪽=B) Logging: - BLE 로그 + CSV에 otsu threshold 값 기록 - 측정 로그에 method label (A/B) 표시 UI: - 설정 패널: floating overlay로 변경 + MaxVolume 슬라이더 추가 - Onboarding: 1페이지로 단순화 (Pager/indicator 제거) - Registration: 한 페이지에 전체 입력 + Skip 버튼 - Placement: 로딩 인디케이터/스캔 텍스트 삭제 Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
This commit is contained in:
@@ -12,8 +12,12 @@ package com.example.medilightv2android.managers
|
||||
* ⚠ 미세 조정 시 이 파일만 수정하면 전체 파이프라인에 반영됨.
|
||||
* 각 상수의 의미와 영향 범위를 아래 주석 참고.
|
||||
*/
|
||||
enum class DetectionMethod { METHOD_A, METHOD_B }
|
||||
|
||||
object GreenZoneConstants {
|
||||
|
||||
var detectionMethod: DetectionMethod = DetectionMethod.METHOD_B
|
||||
|
||||
// ═══════════════════════════════════════════════════════════
|
||||
// 신호 범위
|
||||
// ═══════════════════════════════════════════════════════════
|
||||
|
||||
@@ -266,130 +266,115 @@ private fun segmentToDistancesMm(
|
||||
return Pair(min(d1, d2), max(d1, d2))
|
||||
}
|
||||
|
||||
// ── Parabolic Cap Fitting ──
|
||||
// ── Ellipse Cap Height Fitting ──
|
||||
// bv_estimation.py _ellipse_cap_heights() 1:1 포팅
|
||||
// 8개 경계점(4ch × ant/post)으로 축 정렬 타원 피팅 → cap 높이
|
||||
|
||||
private data class ParabolicCapResult(
|
||||
val s: DoubleArray,
|
||||
val aCap: DoubleArray,
|
||||
val bEff: Double?,
|
||||
val capMode: String,
|
||||
val outlierIdx: Int?
|
||||
private data class EllipseCapResult(
|
||||
val hBot: Double,
|
||||
val hTop: Double,
|
||||
val botKind: String,
|
||||
val topKind: String
|
||||
)
|
||||
|
||||
/**
|
||||
* S(y) 포물선 피팅 → outlier 보정 + b_eff 추출
|
||||
*/
|
||||
private fun parabolicCap(
|
||||
sS: DoubleArray, yS: DoubleArray, aCapS: DoubleArray,
|
||||
outlierSigma: Double = 1.5, r2Min: Double = 0.5
|
||||
): ParabolicCapResult {
|
||||
val n = sS.size
|
||||
val sWork = sS.copyOf()
|
||||
val aWork = aCapS.copyOf()
|
||||
var outlierIdx: Int? = null
|
||||
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>()
|
||||
|
||||
if (n < 3) {
|
||||
return ParabolicCapResult(sWork, aWork, null, "fallback", null)
|
||||
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))
|
||||
}
|
||||
|
||||
// ── Phase 0: Edge-peak 보정 ──
|
||||
val peakIdx = sWork.indices.maxByOrNull { sWork[it] } ?: 0
|
||||
|
||||
if (peakIdx == 0 && n >= 3) {
|
||||
val slope: Double = if (abs(yS[2] - yS[1]) > 1e-6)
|
||||
(sWork[2] - sWork[1]) / (yS[2] - yS[1]) else 0.0
|
||||
val sExtrap = max(sWork[1] + slope * (yS[0] - yS[1]), 1.0)
|
||||
if (sExtrap < sWork[0]) {
|
||||
sWork[0] = sExtrap
|
||||
val ratio = if (sS[0] > 0) sWork[0] / sS[0] else 1.0
|
||||
aWork[0] = aCapS[0] * sqrt(ratio)
|
||||
outlierIdx = 0
|
||||
}
|
||||
} else if (peakIdx == n - 1 && n >= 3) {
|
||||
val slope: Double = if (abs(yS[n - 2] - yS[n - 3]) > 1e-6)
|
||||
(sWork[n - 2] - sWork[n - 3]) / (yS[n - 2] - yS[n - 3]) else 0.0
|
||||
val sExtrap = max(sWork[n - 2] + slope * (yS[n - 1] - yS[n - 2]), 1.0)
|
||||
if (sExtrap < sWork[n - 1]) {
|
||||
sWork[n - 1] = sExtrap
|
||||
val ratio = if (sS[n - 1] > 0) sWork[n - 1] / sS[n - 1] else 1.0
|
||||
aWork[n - 1] = aCapS[n - 1] * sqrt(ratio)
|
||||
outlierIdx = n - 1
|
||||
}
|
||||
if (ptsX.size < 5) {
|
||||
return EllipseCapResult(aCapBot, aCapTop, "fallback", "fallback")
|
||||
}
|
||||
|
||||
// ── Phase 1: 잔차 기반 outlier 보정 (1.5σ) ──
|
||||
val coeffs1 = polyfit2(yS, sWork)
|
||||
val sFitted1 = DoubleArray(n) { polyval2(coeffs1, yS[it]) }
|
||||
val residuals = DoubleArray(n) { sWork[it] - sFitted1[it] }
|
||||
val absRes = DoubleArray(n) { abs(residuals[it]) }
|
||||
val worst = absRes.indices.maxByOrNull { absRes[it] } ?: 0
|
||||
val resStd = std(residuals)
|
||||
val threshold = if (resStd > 1e-6) outlierSigma * resStd else Double.MAX_VALUE
|
||||
// 축 정렬 타원: 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")
|
||||
|
||||
if (worst != outlierIdx && absRes[worst] > threshold) {
|
||||
val sCorr = max(sFitted1[worst], 1.0)
|
||||
sWork[worst] = sCorr
|
||||
val ratio = if (sS[worst] > 0) sWork[worst] / sS[worst] else 1.0
|
||||
aWork[worst] = aCapS[worst] * sqrt(ratio)
|
||||
outlierIdx = worst
|
||||
val (aa, cc, dd, ee) = sol
|
||||
if (aa <= 0 || cc <= 0) {
|
||||
return EllipseCapResult(aCapBot, aCapTop, "fallback", "fallback")
|
||||
}
|
||||
|
||||
// ── Phase 2: 재피팅 → b_eff ──
|
||||
val coeffs2 = polyfit2(yS, sWork)
|
||||
val c2 = coeffs2.first // highest degree coeff
|
||||
val sFitted2 = DoubleArray(n) { polyval2(coeffs2, yS[it]) }
|
||||
val ssRes = (0 until n).sumOf { (sWork[it] - sFitted2[it]).let { d -> d * d } }
|
||||
val meanS = sWork.sum() / n.toDouble()
|
||||
val ssTot = sWork.sumOf { (it - meanS).let { d -> d * d } }
|
||||
val rSq = if (ssTot > 1e-12) 1.0 - ssRes / ssTot else 0.0
|
||||
|
||||
if (c2 < -1e-3 && rSq > r2Min) {
|
||||
val sPeakFit = coeffs2.third - coeffs2.second * coeffs2.second / (4.0 * c2)
|
||||
if (sPeakFit > 0) {
|
||||
val bEff = sqrt(sPeakFit / abs(c2))
|
||||
return ParabolicCapResult(sWork, aWork, bEff, "parabolic", outlierIdx)
|
||||
}
|
||||
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")
|
||||
}
|
||||
|
||||
return ParabolicCapResult(sWork, aWork, null, "fallback", outlierIdx)
|
||||
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")
|
||||
}
|
||||
|
||||
// ── Cap Volume ──
|
||||
/** 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)
|
||||
|
||||
private data class CapResult(val h: Double, val volume: Double, val kind: String)
|
||||
|
||||
/**
|
||||
* 단일 cap (bottom 또는 top) 높이 + 부피 계산
|
||||
* shape: "sphere" = spherical cap, "cone" = 원뿔
|
||||
*/
|
||||
private fun capVolume(
|
||||
sEdge: Double, aCapEdge: Double, bEff: Double?,
|
||||
sMax: Double, capMode: String, shape: String = "sphere"
|
||||
): CapResult {
|
||||
val h: Double
|
||||
var kind: String
|
||||
|
||||
if (capMode == "parabolic" && bEff != null && sMax > 0) {
|
||||
val sRatio = min(sEdge / sMax, 0.999)
|
||||
var hCalc = bEff * (1.0 - sqrt(1.0 - sRatio))
|
||||
hCalc = min(hCalc, aCapEdge) // hemisphere/cone(h=R) 초과 방지
|
||||
h = hCalc
|
||||
kind = "parabolic"
|
||||
} else {
|
||||
h = aCapEdge
|
||||
kind = "fallback"
|
||||
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)
|
||||
}
|
||||
|
||||
val volume: Double
|
||||
if (shape == "cone") {
|
||||
volume = sEdge * h / 3.0
|
||||
kind += " cone"
|
||||
} else {
|
||||
volume = sEdge * h / 2.0 + Math.PI * h * h * h / 6.0
|
||||
kind += " sphere"
|
||||
private fun solve4x4(A: Array<DoubleArray>, 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]
|
||||
}
|
||||
}
|
||||
|
||||
return CapResult(h, volume, kind)
|
||||
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
|
||||
}
|
||||
|
||||
// ── Core BV Computation ──
|
||||
@@ -488,12 +473,12 @@ fun estimateBladderVolume(
|
||||
var aCapS = DoubleArray(order.size) { aCap[order[it]] }
|
||||
val sortedCh = order.map { validChannels[it] }
|
||||
|
||||
// 7) 포물선 피팅 → outlier 보정 + b_eff
|
||||
val paraResult = parabolicCap(sS = sS, yS = yS, aCapS = aCapS)
|
||||
sS = paraResult.s
|
||||
aCapS = paraResult.aCap
|
||||
val bEff = paraResult.bEff
|
||||
val capMode = paraResult.capMode
|
||||
// 7) 타원 피팅 → cap 높이
|
||||
val ellCap = ellipseCapHeights(
|
||||
dAnt, dPost, validChannels,
|
||||
sensorZMm, degreeDeg,
|
||||
aCapS[0], aCapS[n - 1], yS
|
||||
)
|
||||
|
||||
// 8) Core frustum (traditional: h = dy)
|
||||
val dy = DoubleArray(n - 1) { yS[it + 1] - yS[it] }
|
||||
@@ -502,27 +487,12 @@ fun estimateBladderVolume(
|
||||
}
|
||||
val vCore = vFrustum.sum()
|
||||
|
||||
// 9) Caps
|
||||
val sMax = sS.max()
|
||||
|
||||
// Bottom: 항상 sphere
|
||||
val botCap = capVolume(
|
||||
sEdge = sS[0], aCapEdge = aCapS[0], bEff = bEff, sMax = sMax,
|
||||
capMode = capMode, shape = "sphere")
|
||||
|
||||
// Top: cone (정상) / sphere (CH4 short)
|
||||
val topCap = if (edgeIsShort) {
|
||||
capVolume(
|
||||
sEdge = sS[n - 1], aCapEdge = aCapS[n - 1], bEff = bEff, sMax = sMax,
|
||||
capMode = "fallback", shape = "sphere")
|
||||
} else {
|
||||
capVolume(
|
||||
sEdge = sS[n - 1], aCapEdge = aCapS[n - 1], bEff = bEff, sMax = sMax,
|
||||
capMode = capMode, shape = "cone")
|
||||
}
|
||||
// 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
|
||||
|
||||
// 10) 합산
|
||||
val bvMm3 = vCore + botCap.volume + topCap.volume
|
||||
val bvMm3 = vCore + vBottom + vTop
|
||||
|
||||
return BVResult(
|
||||
volumeMl = bvMm3 / 1000.0,
|
||||
@@ -537,12 +507,12 @@ fun estimateBladderVolume(
|
||||
sortedSMm2 = sS,
|
||||
vFrustumMm3 = vFrustum,
|
||||
vCoreMm3 = vCore,
|
||||
vBottomMm3 = botCap.volume,
|
||||
vTopMm3 = topCap.volume,
|
||||
bottomHMm = botCap.h,
|
||||
topHMm = topCap.h,
|
||||
bottomKind = botCap.kind,
|
||||
topKind = topCap.kind,
|
||||
vBottomMm3 = vBottom,
|
||||
vTopMm3 = vTop,
|
||||
bottomHMm = ellCap.hBot,
|
||||
topHMm = ellCap.hTop,
|
||||
bottomKind = ellCap.botKind,
|
||||
topKind = ellCap.topKind,
|
||||
lrRatio = lrRatio,
|
||||
distancePerSample = distancePerSample,
|
||||
delayOffsetMm = delayOffsetMm
|
||||
|
||||
@@ -311,59 +311,78 @@ class PiezoEchoAnalyzer private constructor() {
|
||||
|
||||
// ── Low Echo Detection ──
|
||||
|
||||
/** 적응형 low-echo 임계값: 신호 통계 기반 */
|
||||
fun adaptiveThreshold(signal: DoubleArray): Double {
|
||||
if (signal.size < 10) return lowEchoAmpDefault
|
||||
val sorted = signal.sorted().toDoubleArray()
|
||||
val q25 = sorted[sorted.size / 4]
|
||||
val median = sorted[sorted.size / 2]
|
||||
val q75 = sorted[3 * sorted.size / 4]
|
||||
val iqr = q75 - q25
|
||||
return max(q25, median - 0.5 * iqr)
|
||||
}
|
||||
|
||||
var useAdaptiveThreshold: Boolean = false
|
||||
var useAdaptiveThreshold: Boolean = true
|
||||
var lastOtsuThreshold: Double = 0.0; private set
|
||||
|
||||
/** 단일 1D 채널 → urine region 탐지 */
|
||||
fun detectLowEcho(raw: DoubleArray, denoised: DoubleArray): LowEchoResult? {
|
||||
val sg = denoised
|
||||
if (sg.size < 10) return null
|
||||
|
||||
if (!useAdaptiveThreshold) {
|
||||
return detectLowEchoCore(sg = sg, threshold = lowEchoAmpDefault)
|
||||
val thr = if (useAdaptiveThreshold) {
|
||||
val limit = min(sg.size, postMaxIdx + 1)
|
||||
otsu1d(sg, limit)
|
||||
} else {
|
||||
lowEchoAmpDefault
|
||||
}
|
||||
|
||||
// Adaptive threshold: percentile + depth 보상 조합
|
||||
val thr = computeAdaptiveThreshold(sg)
|
||||
lastOtsuThreshold = thr
|
||||
return detectLowEchoCore(sg = sg, threshold = thr)
|
||||
}
|
||||
|
||||
/**
|
||||
* Adaptive threshold: 신호 통계 + depth attenuation 보상
|
||||
* 1D Otsu threshold — 원은지 연구원 span_utils.py otsu_1d() 1:1 포팅
|
||||
*
|
||||
* 1) Percentile 기반 베이스라인: 전체 신호의 q25~median 사이에서 결정
|
||||
* 2) 초반 peak(피부 반사) 제외: ringSkip(3) 이후 사용
|
||||
* 3) 고정 threshold와의 가중 평균으로 급격한 변동 방지
|
||||
* 히스토그램에서 between-class variance를 최대화하는 임계값 반환.
|
||||
* 저진폭(urine/조직)과 고진폭(벽 echo) 두 모집단을 분리.
|
||||
*/
|
||||
fun computeAdaptiveThreshold(sg: DoubleArray): Double {
|
||||
val skip = GreenZoneConstants.ringSkip
|
||||
val usable = if (sg.size > skip + 10) sg.sliceArray(skip until sg.size) else sg
|
||||
fun otsu1d(values: DoubleArray, limit: Int = values.size, nBins: Int = 64): Double {
|
||||
val n = min(values.size, limit)
|
||||
if (n == 0) return 0.0
|
||||
|
||||
val sorted = usable.sorted()
|
||||
val n = sorted.size
|
||||
val q25 = sorted[n / 4]
|
||||
val q50 = sorted[n / 2]
|
||||
val q75 = sorted[3 * n / 4]
|
||||
val iqr = q75 - q25
|
||||
var vMin = values[0]; var vMax = values[0]
|
||||
for (i in 1 until n) {
|
||||
if (values[i] < vMin) vMin = values[i]
|
||||
if (values[i] > vMax) vMax = values[i]
|
||||
}
|
||||
if (vMax == vMin) return vMin
|
||||
|
||||
// Percentile 기반: 소변 영역은 보통 하위 25~50%
|
||||
val percThr = max(q25 + iqr * 0.3, q50 - iqr * 0.3)
|
||||
val binWidth = (vMax - vMin) / nBins
|
||||
val hist = IntArray(nBins)
|
||||
for (i in 0 until n) {
|
||||
val bin = ((values[i] - vMin) / binWidth).toInt().coerceIn(0, nBins - 1)
|
||||
hist[bin]++
|
||||
}
|
||||
|
||||
// 고정값과 adaptive의 가중 평균 (급격한 변동 방지)
|
||||
val adaptive = percThr.coerceIn(lowEchoAmpDefault * 0.7, lowEchoAmpDefault * 1.5)
|
||||
val blended = lowEchoAmpDefault * 0.4 + adaptive * 0.6
|
||||
val centers = DoubleArray(nBins) { vMin + (it + 0.5) * binWidth }
|
||||
val total = n.toDouble()
|
||||
|
||||
return blended
|
||||
// cumulative probability & mean
|
||||
val cumP = DoubleArray(nBins)
|
||||
val cumMP = DoubleArray(nBins)
|
||||
cumP[0] = hist[0] / total
|
||||
cumMP[0] = cumP[0] * centers[0]
|
||||
for (i in 1 until nBins) {
|
||||
cumP[i] = cumP[i - 1] + hist[i] / total
|
||||
cumMP[i] = cumMP[i - 1] + (hist[i] / total) * centers[i]
|
||||
}
|
||||
val totalM = cumMP[nBins - 1]
|
||||
|
||||
// between-class variance 최대화
|
||||
var bestSigma = -1.0
|
||||
var bestIdx = 0
|
||||
for (t in 0 until nBins - 1) {
|
||||
val w0 = cumP[t]
|
||||
val w1 = 1.0 - w0
|
||||
if (w0 < 1e-6 || w1 < 1e-6) continue
|
||||
val m0 = cumMP[t] / w0
|
||||
val m1 = (totalM - cumMP[t]) / w1
|
||||
val sigmaB = w0 * w1 * (m0 - m1) * (m0 - m1)
|
||||
if (sigmaB > bestSigma) {
|
||||
bestSigma = sigmaB
|
||||
bestIdx = t
|
||||
}
|
||||
}
|
||||
return centers[bestIdx]
|
||||
}
|
||||
|
||||
/**
|
||||
|
||||
@@ -0,0 +1,304 @@
|
||||
package com.example.medilightv2android.managers
|
||||
|
||||
import kotlin.math.abs
|
||||
import kotlin.math.max
|
||||
import kotlin.math.min
|
||||
import kotlin.math.sqrt
|
||||
|
||||
/**
|
||||
* Method A — TVD → SG → Derivative + Plateau-first wall expansion.
|
||||
* 1:1 port of low_echo_detection_method_a.py (원은지 연구원)
|
||||
*
|
||||
* Pipeline:
|
||||
* 1. Condat(2013) 1D TVD
|
||||
* 2. SG smooth
|
||||
* 3. Composite signal y = sg + α|d1| + β|d2|
|
||||
* 4. Sliding plateau score
|
||||
* 5. Plateau span detection + merge
|
||||
* 6. Prominence-based wall selection from plateau edges
|
||||
*/
|
||||
class PiezoEchoAnalyzerA private constructor() {
|
||||
|
||||
companion object {
|
||||
val shared = PiezoEchoAnalyzerA()
|
||||
}
|
||||
|
||||
// ── Parameters (notebook 실측 튜닝값) ──
|
||||
val tvLambda: Double = 3.0
|
||||
val alpha: Double = 1.0
|
||||
val beta: Double = 1.0
|
||||
val postMaxIdx: Int = 80
|
||||
val plateauMinLen: Int = 3
|
||||
val plateauMergeGap: Int = 5
|
||||
val scoreWin: Int = 7
|
||||
val plateauQ: Double = 0.5
|
||||
val edgeDistDecay: Double = 0.15
|
||||
val peakSearchWin: Int = 20
|
||||
val minWallProminence: Double = 3.0
|
||||
val minWallContrastRatio: Double = 1.2
|
||||
val minUrineLen: Int = 3
|
||||
val minEffectiveDist: Int = 3
|
||||
val valleyStopRise: Double = 5.0
|
||||
|
||||
var lastOtsuThreshold: Double = 0.0; private set
|
||||
|
||||
// ── Public API ──
|
||||
|
||||
fun analyzeChannel(rawADC: List<UShort>, channel: Int = 0): ChannelAnalysisResult {
|
||||
val raw = DoubleArray(rawADC.size) { rawADC[it].toDouble() }
|
||||
if (raw.size < 10) {
|
||||
return ChannelAnalysisResult(channel, null, raw, raw.copyOf())
|
||||
}
|
||||
return try {
|
||||
val sg = PiezoEchoAnalyzer.shared.denoise(raw)
|
||||
val result = detectMethodA(raw, sg)
|
||||
ChannelAnalysisResult(channel, result, raw, sg)
|
||||
} catch (_: Exception) {
|
||||
ChannelAnalysisResult(channel, null, raw, raw.copyOf())
|
||||
}
|
||||
}
|
||||
|
||||
fun detectMethodA(raw: DoubleArray, sg: DoubleArray): LowEchoResult? {
|
||||
if (sg.size < 10) return null
|
||||
|
||||
// 1) Composite signal y = sg + α|d1| + β|d2|
|
||||
val d1 = DoubleArray(sg.size - 1) { sg[it + 1] - sg[it] }
|
||||
val d2 = DoubleArray(sg.size - 2) { sg[it + 2] - 2 * sg[it + 1] + sg[it] }
|
||||
val L = d2.size
|
||||
if (L < 10) return null
|
||||
val y = DoubleArray(L) { sg[it] + alpha * abs(d1[it]) + beta * abs(d2[it]) }
|
||||
|
||||
// 2) Plateau score
|
||||
val platScore = slidingScores1d(sg, scoreWin)
|
||||
val limit = min(platScore.size, postMaxIdx + 1)
|
||||
val scoreSorted = platScore.sliceArray(0 until limit).sorted().toDoubleArray()
|
||||
val scoreThr = scoreSorted[(scoreSorted.size * plateauQ).toInt().coerceIn(0, scoreSorted.size - 1)]
|
||||
|
||||
// 3) Plateau spans + merge
|
||||
var plateauSpans = findPlateauSpans(platScore, limit, scoreThr, plateauMinLen)
|
||||
if (plateauSpans.size > 1) {
|
||||
val merged = mutableListOf(plateauSpans[0])
|
||||
for (i in 1 until plateauSpans.size) {
|
||||
val (s, e) = plateauSpans[i]
|
||||
val (ps, pe) = merged.last()
|
||||
if (s - pe <= plateauMergeGap) {
|
||||
merged[merged.size - 1] = Pair(ps, e)
|
||||
} else {
|
||||
merged.add(Pair(s, e))
|
||||
}
|
||||
}
|
||||
plateauSpans = merged
|
||||
}
|
||||
|
||||
// 4) Plateau-first matching
|
||||
for ((spanStart, spanEnd) in plateauSpans) {
|
||||
val leftLo = max(0, spanStart - peakSearchWin)
|
||||
val rightHi = min(sg.size, spanEnd + peakSearchWin + 1)
|
||||
|
||||
val platMean = if (spanEnd >= spanStart) {
|
||||
var sum = 0.0
|
||||
for (i in spanStart..spanEnd) sum += sg[i]
|
||||
sum / (spanEnd - spanStart + 1)
|
||||
} else continue
|
||||
|
||||
if (platMean <= 0) continue
|
||||
|
||||
// Contrast check
|
||||
var leftMax = 0.0
|
||||
for (i in leftLo..min(spanStart, sg.size - 1)) if (sg[i] > leftMax) leftMax = sg[i]
|
||||
var rightMax = 0.0
|
||||
for (i in spanEnd until rightHi) if (sg[i] > rightMax) rightMax = sg[i]
|
||||
if (leftMax / platMean < minWallContrastRatio || rightMax / platMean < minWallContrastRatio) continue
|
||||
|
||||
// Find peaks
|
||||
val leftPks = findPeaksInRange(sg, leftLo, spanStart).filter { it < spanStart }
|
||||
val rightPks = findPeaksInRange(sg, spanEnd, rightHi - 1).filter { it > spanEnd && it <= postMaxIdx }
|
||||
if (leftPks.isEmpty() || rightPks.isEmpty()) continue
|
||||
|
||||
// Score left candidates
|
||||
var bestLeftScore = -1.0; var bestLeft = -1
|
||||
for (p in leftPks) {
|
||||
val prom = sg[p] - findRightValley(sg, p)
|
||||
if (prom < minWallProminence) continue
|
||||
val dist = spanStart - p
|
||||
val effDist = max(dist, minEffectiveDist)
|
||||
val score = prom / (1.0 + edgeDistDecay * effDist)
|
||||
if (score > bestLeftScore) { bestLeftScore = score; bestLeft = p }
|
||||
}
|
||||
|
||||
// Score right candidates
|
||||
var bestRightScore = -1.0; var bestRight = -1
|
||||
for (p in rightPks) {
|
||||
val prom = sg[p] - findLeftValley(sg, p)
|
||||
if (prom < minWallProminence) continue
|
||||
val dist = p - spanEnd
|
||||
val effDist = max(dist, minEffectiveDist)
|
||||
val score = prom / (1.0 + edgeDistDecay * effDist)
|
||||
if (score > bestRightScore) { bestRightScore = score; bestRight = p }
|
||||
}
|
||||
|
||||
if (bestLeft < 0 || bestRight < 0 || bestRight <= bestLeft) continue
|
||||
val urineLen = bestRight - bestLeft
|
||||
if (urineLen < minUrineLen) continue
|
||||
|
||||
val lowSlice = sg.sliceArray(spanStart..spanEnd)
|
||||
val lowMean = lowSlice.average()
|
||||
val wallScore = (sg[bestLeft] + sg[bestRight]) / 2.0 - lowMean
|
||||
|
||||
return LowEchoResult(
|
||||
ant = bestLeft, post = bestRight,
|
||||
lowStart = spanStart, lowEnd = spanEnd,
|
||||
lowMean = lowMean, urineLen = urineLen,
|
||||
score = wallScore * urineLen,
|
||||
innerPeaks = PiezoEchoAnalyzer.shared.findInnerPeaks(sg, bestLeft, bestRight)
|
||||
)
|
||||
}
|
||||
return null
|
||||
}
|
||||
|
||||
// ── Sliding plateau score ──
|
||||
|
||||
private fun slidingScores1d(x: DoubleArray, win: Int): DoubleArray {
|
||||
val w = if (win % 2 == 0) win + 1 else win
|
||||
val half = w / 2
|
||||
val T = x.size
|
||||
val flat = DoubleArray(T)
|
||||
val slope = DoubleArray(T)
|
||||
val low = DoubleArray(T)
|
||||
|
||||
val tt = DoubleArray(w) { (it - half).toDouble() }
|
||||
var ttSqSum = 0.0
|
||||
for (t in tt) ttSqSum += t * t
|
||||
ttSqSum += 1e-12
|
||||
|
||||
for (i in 0 until T) {
|
||||
val wStart = max(0, i - half)
|
||||
val wEnd = min(T - 1, i + half)
|
||||
val wLen = wEnd - wStart + 1
|
||||
|
||||
var sum = 0.0; var sqSum = 0.0
|
||||
val vals = mutableListOf<Double>()
|
||||
for (j in wStart..wEnd) { sum += x[j]; sqSum += x[j] * x[j]; vals.add(x[j]) }
|
||||
val mean = sum / wLen
|
||||
flat[i] = sqrt(max(0.0, sqSum / wLen - mean * mean))
|
||||
|
||||
vals.sort()
|
||||
low[i] = vals[(wLen * 0.2).toInt().coerceIn(0, wLen - 1)]
|
||||
|
||||
var slopeNum = 0.0
|
||||
for (j in wStart..wEnd) {
|
||||
slopeNum += (j - i).toDouble() * (x[j] - mean)
|
||||
}
|
||||
slope[i] = abs(slopeNum / ttSqSum)
|
||||
}
|
||||
|
||||
return DoubleArray(T) { robustZ(low)[it] + robustZ(flat)[it] + robustZ(slope)[it] }
|
||||
}
|
||||
|
||||
private fun robustZ(a: DoubleArray): DoubleArray {
|
||||
val sorted = a.sorted().toDoubleArray()
|
||||
val med = sorted[sorted.size / 2]
|
||||
val deviations = DoubleArray(a.size) { abs(a[it] - med) }
|
||||
val devSorted = deviations.sorted().toDoubleArray()
|
||||
val mad = devSorted[devSorted.size / 2] + 1e-12
|
||||
return DoubleArray(a.size) { (a[it] - med) / (1.4826 * mad) }
|
||||
}
|
||||
|
||||
// ── Plateau span detection ──
|
||||
|
||||
private fun findPlateauSpans(score: DoubleArray, limit: Int, thr: Double, minLen: Int): MutableList<Pair<Int, Int>> {
|
||||
val spans = mutableListOf<Pair<Int, Int>>()
|
||||
var runStart: Int? = null
|
||||
for (i in 0 until limit) {
|
||||
if (score[i] < thr) {
|
||||
if (runStart == null) runStart = i
|
||||
} else {
|
||||
if (runStart != null && (i - runStart) >= minLen) {
|
||||
spans.add(Pair(runStart, i - 1))
|
||||
}
|
||||
runStart = null
|
||||
}
|
||||
}
|
||||
if (runStart != null && (limit - runStart) >= minLen) {
|
||||
spans.add(Pair(runStart, limit - 1))
|
||||
}
|
||||
return spans
|
||||
}
|
||||
|
||||
// ── Cross-validation (채널간 벽 보정) ──
|
||||
// 1:1 port of _cross_validate_walls()
|
||||
|
||||
private val crossValMaxDev = 5
|
||||
private val crossValMaxGradient = 3
|
||||
|
||||
fun crossValidateWalls(
|
||||
walls: List<Pair<Int, Int>?>,
|
||||
centerCh: List<Int> = listOf(0, 1, 2, 3)
|
||||
): List<Pair<Int, Int>?> {
|
||||
val corrected = walls.toMutableList()
|
||||
|
||||
// 1단계: 글로벌 중앙값 기반
|
||||
val ants = corrected.mapNotNull { it?.first }
|
||||
val posts = corrected.mapNotNull { it?.second }
|
||||
if (ants.size >= 3) {
|
||||
val medAnt = ants.sorted()[ants.size / 2]
|
||||
val medPost = posts.sorted()[posts.size / 2]
|
||||
for (i in corrected.indices) {
|
||||
val w = corrected[i] ?: continue
|
||||
val a = if (abs(w.first - medAnt) >= crossValMaxDev) medAnt else w.first
|
||||
val p = if (abs(w.second - medPost) >= crossValMaxDev) medPost else w.second
|
||||
corrected[i] = Pair(a, p)
|
||||
}
|
||||
}
|
||||
|
||||
// 2단계: 인접 center 채널 smoothness
|
||||
val validCenter = centerCh.filter { it < corrected.size && corrected[it] != null }
|
||||
if (validCenter.size < 3) return corrected
|
||||
|
||||
for (field in listOf("ant", "post")) {
|
||||
val vals = validCenter.map { if (field == "ant") corrected[it]!!.first else corrected[it]!!.second }.toMutableList()
|
||||
for (j in vals.indices) {
|
||||
val neighbors = mutableListOf<Int>()
|
||||
if (j > 0) neighbors.add(vals[j - 1])
|
||||
if (j < vals.size - 1) neighbors.add(vals[j + 1])
|
||||
if (neighbors.isEmpty()) continue
|
||||
val neighborMean = neighbors.average()
|
||||
if (abs(vals[j] - neighborMean) > crossValMaxGradient) {
|
||||
val newVal = neighborMean.toInt()
|
||||
val chIdx = validCenter[j]
|
||||
val w = corrected[chIdx]!!
|
||||
corrected[chIdx] = if (field == "ant") Pair(newVal, w.second) else Pair(w.first, newVal)
|
||||
vals[j] = newVal
|
||||
}
|
||||
}
|
||||
}
|
||||
return corrected
|
||||
}
|
||||
|
||||
// ── Peak / Valley helpers ──
|
||||
|
||||
private fun findPeaksInRange(sg: DoubleArray, from: Int, to: Int): List<Int> {
|
||||
val lo = max(0, from); val hi = min(sg.size - 1, to)
|
||||
if (hi - lo < 2) return emptyList()
|
||||
val seg = sg.sliceArray(lo..hi)
|
||||
return PiezoEchoAnalyzer.shared.findPeaks1D(seg).map { it + lo }
|
||||
}
|
||||
|
||||
private fun findRightValley(sg: DoubleArray, peakIdx: Int, maxDist: Int = 20): Double {
|
||||
var v = sg[peakIdx]
|
||||
for (i in (peakIdx + 1) until min(sg.size, peakIdx + maxDist)) {
|
||||
if (sg[i] < v) v = sg[i]
|
||||
else if (sg[i] > v + valleyStopRise) break
|
||||
}
|
||||
return v
|
||||
}
|
||||
|
||||
private fun findLeftValley(sg: DoubleArray, peakIdx: Int, maxDist: Int = 20): Double {
|
||||
var v = sg[peakIdx]
|
||||
for (i in (peakIdx - 1) downTo max(0, peakIdx - maxDist)) {
|
||||
if (sg[i] < v) v = sg[i]
|
||||
else if (sg[i] > v + valleyStopRise) break
|
||||
}
|
||||
return v
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user