refactor: 데모 브랜치 패키지 통일 com.example.medilightv2android → com.medithings.vesiscan

사용자앱 (feature/tab-navigation) 과 패키지 이름 일치. namespace + applicationId
둘 다 com.medithings.vesiscan 로 변경 (기존 데모 앱은 재설치 필요).

## 변경 범위
- Kotlin 148 파일: package + import 문 (714 occurrences)
- 디렉토리 이동: com/example/medilightv2android → com/medithings/vesiscan
  (main, test, androidTest 각각)
- app/build.gradle.kts: namespace, applicationId
- docs/FLAVOR_DEMO_STABLE.md: 참조 갱신
- V41DetectorCH4Test: BvDispatchResult.methodChosen → method (dto field name fix)

## 주의
- applicationId 가 바뀌므로 기존 데모 앱 (com.example.medilightv2android.demo) 은
  Android 관점에서 다른 앱으로 취급 — 재설치 시 PIN/설정 초기화됨.
- Fresh install 권장. 기존 앱 (com.example...) 은 별도로 uninstall 필요.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
This commit is contained in:
2026-07-02 14:30:43 +09:00
parent 55d55225ab
commit bcbc7b440c
159 changed files with 2865 additions and 717 deletions
@@ -0,0 +1,92 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Stability check on V41 detection outputs. Complements SweepStabilizer:
* the stabilizer guards the INPUT, this guards the OUTPUT.
*
* Per-channel rules (relative to the previous live frame):
* • antIdx jump > maxIndexJump → unstable
* • postIdx jump > maxIndexJump → unstable
* • |chord change| / chord > maxChordRatio → unstable
*
* "unstable" is a flag for the diagnostic log, NOT a hard reject — V41's
* own anatomical gate already drops genuinely-bad detections. This layer
* just surfaces "this channel's reading just jumped, double-check".
*/
package com.medithings.vesiscan.walldetect
import com.medithings.vesiscan.walldetect.dto.ChannelResult
import kotlin.math.abs
class DetectionSanity(
private val maxIndexJump: Int = 8,
private val maxChordRatio: Double = 0.20,
private val channels: Int = 6,
) {
private val lastAnt = IntArray(channels) { -1 }
private val lastPost = IntArray(channels) { -1 }
data class ChannelReport(
val channel: Int,
val antJump: Int?,
val postJump: Int?,
val chordChangeRatio: Double?,
val unstable: Boolean,
)
fun reset() {
for (i in lastAnt.indices) { lastAnt[i] = -1; lastPost[i] = -1 }
}
fun check(perChannel: List<ChannelResult>): List<ChannelReport> {
val reports = mutableListOf<ChannelReport>()
for ((ch, cr) in perChannel.withIndex()) {
val ant = cr.antIdx
val post = cr.postIdx
val prevA = lastAnt.getOrElse(ch) { -1 }
val prevP = lastPost.getOrElse(ch) { -1 }
var antJump: Int? = null
var postJump: Int? = null
var chordRatio: Double? = null
var unstable = false
if (ant != null && prevA >= 0) {
val jump = abs(ant - prevA)
antJump = jump
if (jump > maxIndexJump) unstable = true
}
if (post != null && prevP >= 0) {
val jump = abs(post - prevP)
postJump = jump
if (jump > maxIndexJump) unstable = true
}
if (ant != null && post != null && prevA >= 0 && prevP >= 0) {
val curChord = (post - ant).toDouble()
val prevChord = (prevP - prevA).toDouble()
if (prevChord > 0.0) {
val r = abs(curChord - prevChord) / prevChord
chordRatio = r
if (r > maxChordRatio) unstable = true
}
}
reports += ChannelReport(
channel = ch,
antJump = antJump,
postJump = postJump,
chordChangeRatio = chordRatio,
unstable = unstable,
)
// Update state ONLY when we have a valid current detection,
// otherwise comparison against -1 sentinel resumes after gap.
if (ant != null) lastAnt[ch] = ant
if (post != null) lastPost[ch] = post
}
return reports
}
}
@@ -0,0 +1,12 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
* All rights reserved.
*/
package com.medithings.vesiscan.walldetect
object DetectorIds {
const val V2 = "v2"
const val V4_1 = "v4_1"
}
@@ -0,0 +1,219 @@
/*
* Method D detector — port of method_d/detector.py.
*
* 단일 채널 raw → ant/post wall index + subsample refine + diagnostics.
*
* Stage 1 : SG heavy / light denoise (호출부에서 미리 줘도 됨)
* Stage 2 : span (heavy primary, light fallback)
* Stage 3 : ant/post wall peak (find_wall_peak_local)
* Stage 3.5: wall-lumen ratio gate + 최대 2회 recovery (outward stronger peak 탐색)
* Stage 4 : TGC-FP gate (raw post ratio ≥ min_post_raw_ratio)
* Stage 5 : subsample refine (parabolic for peak / d2 fit for shoulder)
*
* 인터페이스는 WallDetector 와 무관 — alignment.py 가 channels[i].urine_len 만
* 필요로 하므로 가벼운 단일-함수 API 로 유지.
*/
package com.medithings.vesiscan.walldetect
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDParams
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDPreprocessing
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDSpan
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDWallSelect
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDWallSelect.CType
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDWallSelect.Side
import kotlin.math.max
data class MethodDResult(
// ant / antRefined / urineLen 는 02bed02 의 cross-channel _apply_ant_tiebreak 가
// 사후에 갱신하므로 var. Python `r.ant = pick[0]`, `r.ant_refined = float(pick[0])`,
// `r.urine_len = int(r.post - pick[0] - 1)` 와 1:1 매칭.
var ant: Int,
// 2026-07-02: post, postRefined, postType 도 var — cross-channel `apply_neighbor_top_validate`
// 와 `apply_inward_post` (piezophantomtest 4830e7c / 5576a59) 가 사후 수정.
var post: Int,
var antRefined: Double,
var postRefined: Double,
val lowStart: Int,
val lowEnd: Int,
val lowAmp: Double,
val inwardWalkAnt: Int,
val inwardWalkPost: Int,
var urineLen: Int,
val antProm: Double,
val postProm: Double,
val antType: CType,
var postType: CType,
val sgHeavy: DoubleArray,
val sgLight: DoubleArray,
/** 초기 ant 후보 [(idx, score), ...] score 내림차순 — 02bed02 cross-channel tie-break 용.
* recovery 전 wall_select 결과를 그대로 노출 (Python: detector.py `ant_candidates`). */
val antCandidates: List<Pair<Int, Double>> = emptyList(),
)
object MethodDDetector {
private fun wallGateOk(antAmp: Double, postAmp: Double, lumenMin: Double, p: MethodDParams): Boolean {
val wallAmp = if (p.wallRatioUsePostOnly) postAmp else kotlin.math.min(antAmp, postAmp)
return wallAmp / max(lumenMin, 1.0) >= p.minWallLumenRatio
}
fun detect(
raw: DoubleArray,
denoisedHeavy: DoubleArray? = null,
denoisedLight: DoubleArray? = null,
otsuRatioOverride: Double? = null,
params: MethodDParams = MethodDParams.DEFAULT,
): MethodDResult? {
val sgHeavy = denoisedHeavy ?: MethodDPreprocessing.preprocessHeavy(raw, params)
val sgLight = denoisedLight ?: MethodDPreprocessing.preprocessLight(raw)
val ratio = otsuRatioOverride ?: params.otsuRatio
// Stage 2: span
val spanRes = MethodDSpan.extractLowEchoSpanWithFallback(sgHeavy, sgLight, ratio, params)
?: return null
val lowStart = spanRes.lowStart
val lowEnd = spanRes.lowEnd
val lowAmp = spanRes.lowAmp
val inwardAnt = spanRes.inwardWalkAnt
val inwardPost = spanRes.inwardWalkPost
// lumen_min — Stage 3.5 wall_lumen_ratio gate 에서 사용
var lumenMin = Double.POSITIVE_INFINITY
for (i in lowStart..lowEnd) if (sgHeavy[i] < lumenMin) lumenMin = sgHeavy[i]
if (!lumenMin.isFinite()) lumenMin = 1.0
// Stage 3: ant + post peak
val antRes = MethodDWallSelect.findWallPeakLocal(
sgLight, lowStart, lowEnd, Side.ANT, params,
dMaxOverride = params.antDMax, inwardWalk = inwardAnt,
)
val postRes = MethodDWallSelect.findWallPeakLocal(
sgLight, lowStart, lowEnd, Side.POST, params,
dMaxOverride = params.dMax, inwardWalk = inwardPost,
)
if (antRes.best == null || postRes.best == null) return null
// 02bed02: 초기 ant 후보 (recovery 전, score 내림차순) — cross-channel tie-break 용.
// Python: tuple((int(c[0]), float(c[4])) for c in sorted(ant_res['candidates'], key=-x[4]))
val antCandidates: List<Pair<Int, Double>> =
antRes.candidates.sortedByDescending { it.score }.map { it.idx to it.score }
var antIdx = antRes.best.idx
var antProm = antRes.best.prom
var antType = antRes.best.type
var postIdx = postRes.best.idx
var postProm = postRes.best.prom
var postType = postRes.best.type
if (postIdx <= antIdx) return null
var urineLen = postIdx - antIdx - 1
if (urineLen < params.minUrineLen) return null
var antAmp = sgLight[antIdx]
var postAmp = sgLight[postIdx]
// Stage 3.5: recovery (최대 2회, 각 side 1회씩)
repeat(2) {
if (wallGateOk(antAmp, postAmp, lumenMin, params)) return@repeat
val side: Side
val curIdx: Int
val curAmp: Double
val otherAmp: Double
val walkKw: Int
val extDMax: Int
if (postAmp <= antAmp) {
side = Side.POST
curIdx = postIdx
curAmp = postAmp
otherAmp = antAmp
walkKw = inwardPost
extDMax = if (params.recoveryExtendOutward)
max(params.dMax, params.postMaxIdx - lowEnd) else params.dMax
} else {
side = Side.ANT
curIdx = antIdx
curAmp = antAmp
otherAmp = postAmp
walkKw = inwardAnt
extDMax = if (params.recoveryExtendOutward)
max(params.antDMax, lowStart) else params.antDMax
}
val ext = MethodDWallSelect.findWallPeakLocal(
sgLight, lowStart, lowEnd, side, params,
dMaxOverride = extDMax, inwardWalk = walkKw,
)
// outward stronger peak 만 후보. closest 우선 (post=오름차순/ant=내림차순).
val sorted = if (side == Side.POST) ext.candidates.sortedBy { it.idx }
else ext.candidates.sortedByDescending { it.idx }
var passing: MethodDWallSelect.Candidate? = null // ratio 통과시키는 closest
var fallback: MethodDWallSelect.Candidate? = null // 통과 못해도 stronger 한 첫 후보
for (c in sorted) {
val outward = if (side == Side.POST) c.idx > curIdx else c.idx < curIdx
if (!outward) continue
val cAmp = sgLight[c.idx]
if (cAmp <= curAmp) continue
if (fallback == null) fallback = c
val minAmp = kotlin.math.min(cAmp, otherAmp)
if (minAmp / max(lumenMin, 1.0) >= params.minWallLumenRatio) {
passing = c
break
}
}
val recovered = passing ?: fallback ?: return null
if (side == Side.POST) {
postIdx = recovered.idx
postProm = recovered.prom
postType = recovered.type
postAmp = sgLight[postIdx]
} else {
antIdx = recovered.idx
antProm = recovered.prom
antType = recovered.type
antAmp = sgLight[antIdx]
}
}
if (!wallGateOk(antAmp, postAmp, lumenMin, params)) return null
// Stage 4: TGC-FP gate — raw post / raw lumen_min ≥ min_post_raw_ratio
val rawHeavy = MethodDPreprocessing.preprocessHeavy(raw, params)
val rawLight = MethodDPreprocessing.preprocessLight(raw)
var rawLumenMin = Double.POSITIVE_INFINITY
for (i in lowStart..lowEnd) if (rawHeavy[i] < rawLumenMin) rawLumenMin = rawHeavy[i]
if (!rawLumenMin.isFinite()) rawLumenMin = 1.0
val rawPostRatio = rawLight[postIdx] / max(rawLumenMin, 1.0)
if (rawPostRatio < params.minPostRawRatio) return null
urineLen = postIdx - antIdx - 1
if (urineLen < params.minUrineLen) return null
// Stage 5: subsample refine
val antRefined = if (antType == CType.PEAK)
MethodDWallSelect.refineParabolic(sgLight, antIdx, true)
else MethodDWallSelect.refineShoulder(sgLight, antIdx)
val postRefined = if (postType == CType.PEAK)
MethodDWallSelect.refineParabolic(sgLight, postIdx, true)
else MethodDWallSelect.refineShoulder(sgLight, postIdx)
return MethodDResult(
ant = antIdx,
post = postIdx,
antRefined = antRefined,
postRefined = postRefined,
lowStart = lowStart,
lowEnd = lowEnd,
lowAmp = lowAmp,
inwardWalkAnt = inwardAnt,
inwardWalkPost = inwardPost,
urineLen = urineLen,
antProm = antProm,
postProm = postProm,
antType = antType,
postType = postType,
sgHeavy = sgHeavy,
sgLight = sgLight,
antCandidates = antCandidates,
)
}
}
@@ -0,0 +1,434 @@
/*
* Method D multichannel runner + cross-channel corrections.
*
* 1:1 port of piezophantomtest `library/runners.py::method_d()` +
* `library/cross_channel.py` (2026-07-02 4830e7c + 5576a59)
*
* 파이프라인:
* 1) heavy = SG(7,3) + oscfar_median(win=9, iter=4) ← 02bed02 win 5→9
* light = SG(7,3)
* 2) apply_tgc_pipeline(heavy / light) — center_ch=None → all channels
* 3) 채널별 otsu_ratio × cos(beam_angle) ← PiezoHW.degreeAll
* 4) MethodDDetector.detect()
* 5) cross-channel 후처리 (순서 고정):
* ① applyAntTiebreak — 전벽 교차보정 (4830e7c 신 버전)
* ② applyNeighborTopValidate — 최상단 채널 FP/FN 검증 (신규 4830e7c)
* ③ applyInwardPost — 최상단 widest post 과확장 교정 (신규 5576a59)
*
* 출력 컨트랙트: alignment.py 가 `dets[i].urine_len` 만 의존 →
* MethodDResult.urineLen 또는 null 의 List 로 충분.
*/
package com.medithings.vesiscan.walldetect
import com.medithings.vesiscan.managers.AlignmentConstants
import com.medithings.vesiscan.managers.PiezoHW
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDParams
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDPreprocessing
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDTgc
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDWallSelect
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDWallSelect.CType
import com.medithings.vesiscan.walldetect.algo.methodd.MethodDWallSelect.Side
import kotlin.math.abs
import kotlin.math.cos
import kotlin.math.min as kmin
object MethodDRunner {
// ── cross-channel tuning constants (cross_channel.py:30-45) ────────────────
/** 깊은쪽 교정 시 "합의에 이 이상 더 가까운 후보" 요구 (곡률 보존). */
private const val ANT_OUTLIER_IMPROVE_MARGIN_MM = 5.0
/** 교정 대상이 합의에서 이 이상 떨어지면 보정 보류 (얕은쪽은 제외). */
private const val ANT_CONSENSUS_CLOSE_TOL_MM = 6.0
/** dominant-pick 보호: raw 픽 score >= 이 배수 × 차순위 → 합의 교정에서 보존. */
private const val ANT_DOMINANT_RATIO = 2.0
/** neighbor_top_validate: nb-1→nb slope 가 이 이하면 "비증가 추세". */
private const val NTV_SLOPE_EPS = 2.0
/** neighbor_top_validate: top 이 예측보다 이만큼 깊으면 위반 (FP). */
private const val NTV_DEV_DELTA = 6.0
/** cross-channel 3개 스위치 (Python 과 동일 default). */
private const val NEIGHBOR_TOP_VALIDATE = true
private const val INWARD_POST_SEARCH = true
/**
* 6채널 (또는 N채널) raw 신호 → 채널별 MethodDResult? 리스트.
* raw[ch] 길이가 다르면 그대로 처리 (각 채널 독립).
*/
fun detectMultichannel(
signals: List<DoubleArray>,
params: MethodDParams = MethodDParams.DEFAULT,
beamAnglesDeg: DoubleArray? = null,
applyTgc: Boolean = true,
): List<MethodDResult?> {
val angles = beamAnglesDeg ?: PiezoHW.degreeAll
// 1) per-channel SG denoise (heavy + light)
val heavyList = signals.map { MethodDPreprocessing.preprocessHeavy(it, params) }
val lightList = signals.map { MethodDPreprocessing.preprocessLight(it) }
// 2) TGC per channel (Python apply_tgc_pipeline default center_ch=None → all)
val heavyTgc = if (applyTgc) MethodDTgc.applyTgcPipeline(heavyList) else heavyList
val lightTgc = if (applyTgc) MethodDTgc.applyTgcPipeline(lightList) else lightList
// 3+4) per-channel cos-angle adjusted otsu_ratio + detect
val results = MutableList(signals.size) { ch ->
val angleDeg = if (ch < angles.size) angles[ch] else 0.0
val chRatio = params.otsuRatio * cos(Math.toRadians(angleDeg))
MethodDDetector.detect(
raw = signals[ch],
denoisedHeavy = heavyTgc[ch],
denoisedLight = lightTgc[ch],
otsuRatioOverride = chRatio,
params = params,
)
}
// 5) cross-channel 후처리 (Python `apply_cross_channel` 순서 고정)
applyCrossChannel(results, angles, lightTgc, heavyTgc, signals, params)
return results
}
/**
* Cross-channel 3-stage 후처리 (Python `apply_cross_channel`).
* ① ant_tiebreak → ② neighbor_top_validate → ③ inward_post
*/
private fun applyCrossChannel(
results: MutableList<MethodDResult?>,
angles: DoubleArray,
light: List<DoubleArray>,
heavy: List<DoubleArray>,
raw6ch: List<DoubleArray>,
params: MethodDParams,
) {
applyAntTiebreak(results, angles)
if (NEIGHBOR_TOP_VALIDATE) {
applyNeighborTopValidate(results, angles, light, heavy, raw6ch, params)
}
if (INWARD_POST_SEARCH) {
applyInwardPost(results, angles)
}
}
// ─────────────────────────────────────────────────────────────────────────
// ① ant_tiebreak — 전벽 교차보정 (2026-07-02: 4830e7c 신 버전)
// ─────────────────────────────────────────────────────────────────────────
/**
* 이웃 center 채널 전벽 z 합의 (LOO median) 기준 전벽 보정 (in-place).
*
* 단일 임계 `devTol=8mm` + 방향 분기 (구 nf_tol 8 / outlier_tol 9 통합):
* [최우선] dominant-pick 보호: raw 픽 score >= 2.0 × 차순위 → 보존
* (1) 동률(top1/top2 < tieRatio): 동률 band 중 합의 최근접 선택
* (2) 얕은쪽 (합의 - dev_tol 미만 = near-field 아티팩트): 무조건 재선택
* (3) 깊은쪽 (합의 + dev_tol 초과): improve_margin 이상 더 가까운 후보 있을 때만
*
* 얕은쪽 제외 후 합의 최근접 pick. 합의-근접 가드 (얕은쪽 교정은 항상, 그 외엔 pick 이
* 합의 근접일 때만).
*/
private fun applyAntTiebreak(
results: MutableList<MethodDResult?>,
angles: DoubleArray,
tieRatio: Double = AlignmentConstants.ANT_TIEBREAK_RATIO,
devTol: Double = AlignmentConstants.ANT_NEARFIELD_TOL_MM,
improveMargin: Double = ANT_OUTLIER_IMPROVE_MARGIN_MM,
) {
if (results.isEmpty()) return
// 현재 center 채널 ant z 좌표 (LOO median 의 base)
val baseZ = HashMap<Int, Double>()
for (i in AlignmentConstants.CENTER_CH) {
if (i < results.size) {
val r = results[i] ?: continue
baseZ[i] = zOf(angles, i, r.ant)
}
}
for (i in AlignmentConstants.CENTER_CH) {
if (i >= results.size) continue
val r = results[i] ?: continue
val cands = r.antCandidates
if (cands.size < 2) continue
val cons = AlignmentConstants.CENTER_CH
.filter { it != i }
.mapNotNull { baseZ[it] }
if (cons.isEmpty()) continue
val cz = median(cons)
// dominant-pick 보호 (최우선)
if (cands[0].second >= ANT_DOMINANT_RATIO * cands[1].second
&& r.ant == cands[0].first) {
continue
}
val cur = zOf(angles, i, r.ant)
val curDev = abs(cur - cz)
val bestDev = cands.minOf { abs(zOf(angles, i, it.first) - cz) }
val isTie = cands[0].second / maxOf(cands[1].second, 1e-9) < tieRatio
val isShallow = cur < cz - devTol
val isDeep = cur > cz + devTol
// 진입 + 후보 pool (방향 분기)
val pool: List<Pair<Int, Double>> = when {
isTie -> {
val thr = cands[0].second / tieRatio
cands.filter { it.second >= thr }
}
isShallow -> cands.toList() // 얕은 아티팩트 → 무조건
isDeep && bestDev < curDev - improveMargin -> cands.toList()
else -> continue // 압승 & 정상 범위 → 유지
}
// 얕은(아티팩트) 후보 제외, 합의 최근접 선택
val eligPrelim = pool.filter { zOf(angles, i, it.first) >= cz - devTol }
val elig = if (eligPrelim.isEmpty()) pool else eligPrelim
val pick = elig.minBy { abs(zOf(angles, i, it.first) - cz) }
// 합의-근접 가드: 얕은쪽 교정은 항상, 그 외엔 pick 이 합의 근접일 때만
val pickDev = abs(zOf(angles, i, pick.first) - cz)
if (pickDev > ANT_CONSENSUS_CLOSE_TOL_MM && !isShallow) continue
if (pick.first != r.ant && r.post > pick.first) {
r.ant = pick.first
r.antRefined = pick.first.toDouble()
r.urineLen = r.post - pick.first - 1
}
}
}
// ─────────────────────────────────────────────────────────────────────────
// ② neighbor_top_validate — 최상단 채널 FP/FN 검증 (신규 4830e7c)
// ─────────────────────────────────────────────────────────────────────────
/**
* 이웃 후위벽 trend 로 최상단 채널 검증 (in-place).
*
* Rule B (trend-FP): 검출 top 후위벽 z 를 아래 두 채널 기울기로 예측. 비증가추세인데
* top 이 예측보다 깊게(>6mm) jump → nb urine span 시드 재탐색 → drop/교체.
* Rule A (FN): 미검출 top 을 바로 아래 검출 채널 span 시드로 재탐색 복원(gate 통과 시).
*/
private fun applyNeighborTopValidate(
results: MutableList<MethodDResult?>,
angles: DoubleArray,
light: List<DoubleArray>,
heavy: List<DoubleArray>,
raw6ch: List<DoubleArray>,
params: MethodDParams,
) {
val n = kmin(4, results.size)
val origNone = (0 until n).filter { results[it] == null }.toHashSet()
val det = (0 until n).filter { results[it] != null }.sorted()
// ---- Rule B ----
if (det.size >= 3) {
val top = det[0]; val nb1 = det[1]; val nb2 = det[2]
val rTop = results[top]!!
val rNb1 = results[nb1]!!
val rNb2 = results[nb2]!!
val zTop = zOf(angles, top, rTop.post)
val zNb1 = zOf(angles, nb1, rNb1.post)
val zNb2 = zOf(angles, nb2, rNb2.post)
val slope = zNb1 - zNb2
val pred = zNb1 + slope
if (slope <= NTV_SLOPE_EPS && (zTop - pred) > NTV_DEV_DELTA) {
val rr = researchPostInSpan(
light[top], heavy[top], raw6ch[top], rNb1.lowStart, rNb1.lowEnd, params
)
if (rr == null || (zOf(angles, top, rr.post) - pred) > NTV_DEV_DELTA) {
results[top] = null
} else {
rTop.post = rr.post
rTop.postRefined = rr.post.toDouble()
rTop.postType = rr.postType
rTop.urineLen = rTop.post - rTop.ant - 1
}
}
}
// ---- Rule A ----
val first = (0 until n).firstOrNull { results[it] != null }
if (first != null && first >= 1 && (first - 1) in origNone) {
val top = first - 1
val nb = first
val rNb = results[nb]!!
val ls = rNb.lowStart
val le = rNb.lowEnd
val rr = researchPostInSpan(light[top], heavy[top], raw6ch[top], ls, le, params)
if (rr != null) {
val lowAmp = minInRange(heavy[top], ls, le)
results[top] = MethodDResult(
ant = rr.ant,
post = rr.post,
antRefined = rr.ant.toDouble(),
postRefined = rr.post.toDouble(),
lowStart = ls,
lowEnd = le,
lowAmp = lowAmp,
inwardWalkAnt = 0,
inwardWalkPost = 0,
urineLen = rr.post - rr.ant - 1,
antProm = 0.0,
postProm = 0.0,
antType = rr.antType,
postType = rr.postType,
sgHeavy = heavy[top],
sgLight = light[top],
antCandidates = emptyList(),
)
}
}
}
/** _research_post_in_span 결과. */
private data class ResearchResult(val ant: Int, val post: Int, val antType: CType, val postType: CType)
/**
* seed span(ls, le) 에서 ant / post 재탐색 + wall/raw gate. 실패 시 null.
* detector 와 동일한 gate (_wall_gate_ok + min_post_raw_ratio) 로 재검증.
*/
private fun researchPostInSpan(
lt: DoubleArray, hv: DoubleArray, raw: DoubleArray,
ls: Int, le: Int, params: MethodDParams,
): ResearchResult? {
if (le <= ls) return null
val antRes = MethodDWallSelect.findWallPeakLocal(
lt, ls, le, Side.ANT, params,
dMaxOverride = params.antDMax, inwardWalk = 0,
)
val postRes = MethodDWallSelect.findWallPeakLocal(
lt, ls, le, Side.POST, params,
dMaxOverride = params.dMax, inwardWalk = 0,
)
val antBest = antRes.best ?: return null
val postBest = postRes.best ?: return null
val ai = antBest.idx
val pi = postBest.idx
if (pi <= ai || (pi - ai - 1) < params.minUrineLen) return null
// wall gate
val lumenMin = minInRange(hv, ls, le)
val wallAmp = if (params.wallRatioUsePostOnly) lt[pi]
else kmin(lt[ai], lt[pi])
if (wallAmp / maxOf(lumenMin, 1.0) < params.minWallLumenRatio) return null
// raw ratio gate (SG-only 재계산)
val rawHv = MethodDPreprocessing.preprocessHeavy(raw, params)
val rawLt = MethodDPreprocessing.preprocessLight(raw)
val rawLumen = minInRange(rawHv, ls, le)
if (rawLt[pi] / maxOf(rawLumen, 1.0) < params.minPostRawRatio) return null
return ResearchResult(ai, pi, antBest.type, postBest.type)
}
// ─────────────────────────────────────────────────────────────────────────
// ③ inward_post — 최상단 widest 후벽 과확장 교정 (신규 5576a59)
// ─────────────────────────────────────────────────────────────────────────
/**
* 최상단 검출 center 채널이 widest (적도가 fan 위) 면, 그 채널 후벽을
* [post-win..post] 범위의 더 inward 한 후보(urine-wall ratio ≥ gate)로 교체.
* (a) 분리형 peak, (b) shoulder.
*/
private fun applyInwardPost(
results: MutableList<MethodDResult?>,
angles: DoubleArray,
win: Int = 16,
ratioGate: Double = 1.15,
frac: Double = 0.75,
) {
val det = (0 until kmin(4, results.size))
.filter { results[it] != null }
.sorted()
if (det.size < 2) return
fun dOf(i: Int): Double {
val r = results[i]!!
val ang = if (i < angles.size) angles[i] else 0.0
return (r.post - r.ant) * cos(Math.toRadians(ang))
}
val s = det.associateWith { dOf(it) * dOf(it) }
val top = det[0]
val topS = s[top]!!
val maxS = s.values.max()
val secondS = s[det[1]]!!
if (topS < maxS || topS < secondS) return
val r = results[top]!!
val p = r.post
val ls = r.lowStart
val le = r.lowEnd
val lt = r.sgLight
val hv = r.sgHeavy
if (le <= ls || p - 2 <= ls + 3) return
val base = minInRange(hv, ls, le)
if (base <= 0 || base.isNaN() || base.isInfinite()) return
val postRatio = if (p in lt.indices) lt[p] / base else 0.0
val lo = maxOf(ls + 3, p - win)
val cands = mutableListOf<Int>()
for (i in lo until p - 2) {
if (i !in lt.indices) continue
if (lt[i] / base < ratioGate) continue
// (a) 분리형 peak: light local max, i~post 사이 valley, 깊은 peak 강도 frac 이상
val isPeak = i - 1 in lt.indices && i + 1 in lt.indices
&& lt[i] >= lt[i - 1] && lt[i] > lt[i + 1]
if (isPeak) {
var minLtIP = lt[i]
for (j in i..p) if (j in lt.indices && lt[j] < minLtIP) minLtIP = lt[j]
if (minLtIP < lt[i] * 0.97 && lt[i] >= postRatio * base * frac) {
cands.add(i)
continue
}
}
// (b) shoulder: 3-샘플 plateau + 앞쪽 상승 + 뒤에 더 깊은 peak
if (i + 3 <= p && i - 2 in lt.indices) {
var maxWin = lt[i]; var minWin = lt[i]
for (j in i until i + 3) {
if (j in lt.indices) {
if (lt[j] > maxWin) maxWin = lt[j]
if (lt[j] < minWin) minWin = lt[j]
}
}
val flat = (maxWin - minWin) < 0.02 * lt[i]
val risingBefore = lt[i] > lt[i - 2] + 0.03 * base
var maxAfter = lt[i + 3]
for (j in i + 3..p) if (j in lt.indices && lt[j] > maxAfter) maxAfter = lt[j]
val higherAfter = maxAfter > lt[i] * 1.03
if (flat && risingBefore && higherAfter) {
cands.add(i)
}
}
}
if (cands.isEmpty()) return
val newPost = cands.min()
if (newPost < p - 3) {
r.post = newPost
r.postRefined = newPost.toDouble()
r.urineLen = r.post - r.ant - 1
}
}
// ─────────────────────────────────────────────────────────────────────────
// 공통 helpers
// ─────────────────────────────────────────────────────────────────────────
/** z(i, idx) = (DELAY_OFFSET_MM + idx * DISTANCE_PER_SAMPLE) * cos(angle_i) — sample_to_ap_depth. */
private fun zOf(angles: DoubleArray, ch: Int, idx: Int): Double {
val ang = if (ch < angles.size) angles[ch] else 0.0
return (PiezoHW.delayOffsetMm + idx * PiezoHW.distancePerSample) * cos(Math.toRadians(ang))
}
/** numpy.median 동작 매칭 — 짝수 길이면 두 가운데 값의 평균. */
private fun median(xs: List<Double>): Double {
if (xs.isEmpty()) return 0.0
val sorted = xs.sorted()
val n = sorted.size
return if (n % 2 == 1) sorted[n / 2]
else (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0
}
/** arr[from..to] 최소 (numpy min 매칭). 유효 범위 밖은 skip. 없으면 0.0. */
private fun minInRange(arr: DoubleArray, from: Int, to: Int): Double {
var m = Double.POSITIVE_INFINITY
val lo = maxOf(0, from)
val hi = kmin(arr.size - 1, to)
for (i in lo..hi) if (arr[i] < m) m = arr[i]
return if (m.isFinite()) m else 0.0
}
}
@@ -0,0 +1,137 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Per-channel sweep stabilization layer.
*
* ── Why ─────────────────────────────────────────────────────
* Firmware is observed to occasionally emit ONE channel's ADC
* stream at ~88-90 % of the normal amplitude (TGC gain register
* glitch / Vref sag / ADC trigger miss; observed 25 % of frames
* in the 2026-04-30 11:47 capture). Each anomaly mis-locates
* ant/post for that channel and corrupts the sphere fit.
*
* ── How ─────────────────────────────────────────────────────
* 1) Maintain a ring buffer of the last `historySize` raw ADC
* sweeps (per channel × 100 samples).
* 2) For each new sweep, per channel:
* head_mean(now) / median(head_mean over history) → ratio
* if |ratio − 1| > tolerance ⇒ ANOMALY
* replace channel samples with element-wise median
* across the history → V41 sees a robust value
* else ⇒ pass through
* 3) Expose the anomalous channel list so the caller (VM) can
* surface it in the mbb diagnostic log.
*
* ── What this does NOT do ───────────────────────────────────
* • Does NOT re-scale the bad channel (no artificial correction).
* • Does NOT touch detection logic — V41 stays bit-identical.
* • Does NOT cross-talk between channels.
*
* The contract: "I either pass the live sample through unchanged,
* or replace it with the temporal median of recent good values."
*/
package com.medithings.vesiscan.walldetect
import kotlin.math.abs
class SweepStabilizer(
private val historySize: Int = 3,
/** Allowed band for current head-mean / running median. Outside → anomaly. */
private val tolerance: Double = 0.08, // ±8 %
/** Number of leading samples used to estimate per-channel "amplitude". */
private val headSize: Int = 8,
private val channels: Int = 6,
) {
/** Snapshot of one full sweep's ADC matrix (channel × samples). */
private val history = ArrayDeque<Array<IntArray>>()
data class Anomaly(
val channel: Int,
/** current head-mean ÷ running-median head-mean. ~1.0 is normal. */
val ratio: Double,
/** true if temporal median was substituted; false if first-frame
* (no history yet → cannot replace, passed through unchanged). */
val replaced: Boolean,
)
data class Result(
val adc: List<List<Int>>,
val anomalies: List<Anomaly>,
)
/** Forget all history — call on disconnect / probe re-positioning so the
* stabilizer doesn't compare new captures against stale data. */
fun reset() = history.clear()
fun submit(adc: List<List<Int>>): Result {
require(adc.size >= channels) { "expected $channels channels, got ${adc.size}" }
// Materialise current sweep into IntArrays for efficient median work.
val current = Array(channels) { ch ->
val src = adc[ch]
IntArray(src.size) { src[it] }
}
// No history yet → pass through. Seed the buffer.
if (history.isEmpty()) {
history.addLast(current.deepCopy())
return Result(adc, emptyList())
}
val anomalies = mutableListOf<Anomaly>()
val outChannels = Array(channels) { ch ->
val curHead = headMean(current[ch])
val histHeads = history.map { headMean(it[ch]) }.sorted()
val medHead = histHeads[histHeads.size / 2]
if (medHead <= 0.0) {
current[ch] // degenerate baseline; can't judge
} else {
val ratio = curHead / medHead
if (abs(ratio - 1.0) > tolerance) {
// Anomaly — replace with element-wise median across history
// (NOT including the suspect current frame).
anomalies += Anomaly(channel = ch, ratio = ratio, replaced = true)
elementWiseMedian(history.map { it[ch] }, current[ch].size)
} else {
current[ch]
}
}
}
// Push into history. We push the STABILIZED version so a single
// anomaly can't pollute the median for the next 3 frames.
history.addLast(outChannels.deepCopy())
while (history.size > historySize) history.removeFirst()
return Result(
adc = outChannels.map { it.toList() },
anomalies = anomalies,
)
}
private fun headMean(samples: IntArray): Double {
val n = minOf(headSize, samples.size)
if (n == 0) return 0.0
var sum = 0L
for (i in 0 until n) sum += samples[i]
return sum.toDouble() / n
}
/** Element-wise median across the given snapshots, all assumed length `len`. */
private fun elementWiseMedian(snapshots: List<IntArray>, len: Int): IntArray {
if (snapshots.isEmpty()) return IntArray(len)
val out = IntArray(len)
val tmp = IntArray(snapshots.size)
for (i in 0 until len) {
for (j in snapshots.indices) tmp[j] = snapshots[j][i]
tmp.sort()
out[i] = tmp[tmp.size / 2]
}
return out
}
private fun Array<IntArray>.deepCopy(): Array<IntArray> =
Array(this.size) { this[it].copyOf() }
}
@@ -0,0 +1,106 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* V2Detector v2.2.0 — aligned to py2/for_app_share/low_echo_detection_method_b.
* Replaces both the JS-V2 lumen-first (PR-2) and the plateau-based v2.1.0
* (PR-11) with the actual SSOT mirrored by JS detect_lumen_first.
*
* Pipeline (method_b):
* 1. SG denoise (window=5, polyorder=2)
* 2. Adaptive low-echo threshold via 1-D Otsu on sg[0..POST_MAX_IDX]
* 3. low_mask = sg ≤ T → contiguous spans (≥ LOW_MIN_LEN=3) → merge gaps
* 4. peak_min = max(low_mean + 30, T) — both criteria required
* 5. wall selection by prominence × edge-distance decay (EDGE_DIST_DECAY=0.12)
* with valley-walk stop rise=50
* 6. post > POST_MAX_IDX → back-half retry
* 7. urine_len ≥ 3 required
*/
package com.medithings.vesiscan.walldetect
import com.medithings.vesiscan.walldetect.algo.DetectLumenFirst
import com.medithings.vesiscan.walldetect.algo.Geometry
import com.medithings.vesiscan.walldetect.algo.SubsampleRefine
import com.medithings.vesiscan.walldetect.dto.ChannelResult
import com.medithings.vesiscan.walldetect.dto.DetectionResult
import com.medithings.vesiscan.walldetect.dto.DetectionSummary
import com.medithings.vesiscan.walldetect.dto.SweepInput
import kotlin.math.sqrt
class V2Detector : WallDetector {
override val algorithmId: String = DetectorIds.V2
/** v2.0.0 (JS lumen-first) → v2.1.0 (plateau, scrapped) → v2.2.0 (py2 method_b). */
override val algorithmVersion: String = "v2.2.0"
override fun detect(input: SweepInput): DetectionResult {
val t0 = System.nanoTime()
val perCh: List<ChannelResult> = (0..5).map { ch -> detectChannel(ch, input) }
val chordsMm: List<Float> = perCh.mapNotNull { it.chordMm }
val (chordMmMean, chordMmStd) = if (chordsMm.isEmpty()) null to null
else {
val mu = chordsMm.average()
val v = chordsMm.map { (it - mu) * (it - mu) }.sum() / chordsMm.size
mu.toFloat() to sqrt(v).toFloat()
}
return DetectionResult(
requestId = input.requestId,
timestampMs = input.timestampMs,
algorithm = algorithmId,
algorithmVersion = algorithmVersion,
processingMs = (System.nanoTime() - t0) / 1_000_000.0,
perChannel = perCh,
summary = DetectionSummary(
matchCount = perCh.count { it.antIdx != null && it.postIdx != null },
chordMmMean = chordMmMean,
chordMmStd = chordMmStd,
),
)
}
private fun detectChannel(ch: Int, input: SweepInput): ChannelResult {
val rawList = input.adc[ch]
val rawArr = IntArray(rawList.size) { rawList[it] }
// py2 method_b: Mode.Otsu — adaptive threshold per channel
val det = DetectLumenFirst.detect(rawArr, DetectLumenFirst.Mode.Otsu)
val antRefined = det.ant?.let {
SubsampleRefine.refineParabolic(det.sg, it, SubsampleRefine.Kind.PEAK)
}
val postRefined = det.post?.let {
SubsampleRefine.refineParabolic(det.sg, it, SubsampleRefine.Kind.PEAK)
}
val antMm = antRefined?.let { Geometry.sampleToMm(it).toFloat() }
?: det.ant?.let { Geometry.sampleToMm(it.toDouble()).toFloat() }
val postMm = postRefined?.let { Geometry.sampleToMm(it).toFloat() }
?: det.post?.let { Geometry.sampleToMm(it.toDouble()).toFloat() }
val chordMm = if (antMm != null && postMm != null) postMm - antMm else null
// method_b uses an adaptive scalar T (Otsu); fill the threshold trace
// with that constant so the chart still renders a horizontal reference.
val thrT = det.adaptiveT.toFloat()
val thrTrace = List(det.sg.size) { thrT }
return ChannelResult(
ch = ch,
sg = det.sg.map { it.toFloat() },
threshold = thrTrace,
antIdx = det.ant,
postIdx = det.post,
antRefined = antRefined?.toFloat(),
postRefined = postRefined?.toFloat(),
antMm = antMm,
postMm = postMm,
lumenStart = det.lowStart,
lumenEnd = det.lowEnd,
chordMm = chordMm,
clipping = null,
v41Diag = null,
)
}
}
@@ -0,0 +1,409 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* ═════════════════════════════════════════════════════════════════════════
* V41Detector — Charles V4.1 algorithm (PUBLIC API for team members)
* ═════════════════════════════════════════════════════════════════════════
*
* ▶ USAGE (one-liner)
* ──────────────────
* val result = V41Detector().detect(sweepInput)
*
* ▶ TYPICAL FLOW
* ─────────────
* val mbb = measureCommands.measureFull().getOrThrow() // BLE mbb command
* val sweep = mbb.toSweepInput(probe = SampleSweeps.defaultProbe())
* val result = V41Detector().detect(sweep)
* val bv = result.summary.bvDispatch?.bvMl // mL
* val gated = result.perChannel.count { it.v41Diag?.gated == true }
*
* ▶ INPUT · `dto.SweepInput`
* ──────
* adc : List<List<Int>> // 6 × 100 ADC matrix
* probe : ProbeProfileDto // metadata only — runtime uses WdProbe
* gainDb : Int // V4.1 ignores; V2 fixed-thr scaling
*
* ▶ OUTPUT · `dto.DetectionResult`
* ───────
* algorithm = "v4_1"
* algorithmVersion = "v4.1.1"
* processingMs : Double
* perChannel[6] : ChannelResult
* ├─ sg[100] sg-smoothed envelope
* ├─ threshold[100] OS-CFAR per-sample T
* ├─ antIdx / postIdx wall indices (null when gated out)
* ├─ antRefined / postRefined parabolic sub-sample refine
* ├─ antMm / postMm converted via Geometry.sampleToMm
* ├─ lumenStart / End low-echo span on sg
* ├─ chordMm postMm − antMm
* ├─ clipping ADC saturation flag
* └─ v41Diag ★ V4.1-only diagnostics (see V41Diagnostics.kt)
* summary : DetectionSummary
* ├─ matchCount / chordMmMean / chordMmStd
* ├─ tier / scoreMean / scoreMin / scoreMax / gatedCount
* ├─ v41Sphere Kasa→LM (≥ 4 gated channels)
* └─ bvDispatch ★ multi-method BV (see BvEstimation.kt)
*
* ▶ ALGORITHM PIPELINE (per channel unless marked sweep)
* ─────────────────────────────────────────────────────
* 1.a Denoising.sgSmooth — Savitzky-Golay (5,2)
* 1.b WaveletDenoise.denoise — Sun 2024 db4 3-level adaptive
* 1.c WaveletDenoise.diagnose — energy ratio L1/L2/L3
* 1.d WaveletDenoise.dwt — coefficients (scalogram)
* 2.a ThresholdOsCfar.perSample — Rohling 1983 (k=0.7, scale=1.05)
* 2.b Clipping.isClipped — ADC saturation flag
* 2.c PeakDetection.findPeaks1D — strict local maxima
* 2.d SpanUtils — boolean span merging
* 3.a DetectLumenFirst (Mode.Adaptive) — lumen → wall-pair selection
* 3.b WallSelect.selectWallByProminence— prominence × edge-distance decay
* 3.c SubsampleRefine.refineParabolic — Cespedes 1995
* 4.a ContrastAux.farPostContrast — far-post recovery (s_contrast)
* 4.b BModeScore.score — composite [0,1] (RAW envelope!)
* 4.c AnatomicalGate.apply — PHANTOM_530 [22,60] mm depth
* 6.b BvFromSphere.bvFromChord — single-channel chord-as-D BV
* ─── sweep level ───
* 5.a Geometry.buildWallPointsVisual — visual coords (≤ 12 wall pts)
* 5.d SphereFit2Step.fit2Step — Kasa → LM (auto sphere/circle)
* 6.a BvFromSphere.bvFromSphere — 4/3 π R³ / 1000
* ★ BvEstimation.estimate — multi-method BV dispatch
*
* ▶ KEY CONFIG (`core.WdConfig`)
* ────────────
* DPS_DEFAULT = 1.9309 mm/sample (amode_simulator V4 SSOT)
* DELAY_MM_DEFAULT = 6.85 mm
* LOW_ECHO_AMP = 1250 ADC (V2 fixed; V4.1 uses Otsu)
* PHANTOM_530 = R 50.20 mm, BV 530 mL, ant 32 mm
* WdProbe = v2 (30° device, Snell-refracted angles) by default
*
* ▶ DEPENDENCIES (do not remove)
* ─────────────
* data/protocol/ — packet build/parse + CRC16
* data/command/MeasureCommands.measureFull() → mbb call
* domain/walldetect/algo/ — 17+ algorithm modules (see Pipeline above)
* domain/walldetect/dto/ — kotlinx.serialization DTOs
* data/ble/BleConnector — BLE GATT connection
*
* ▶ NOTES for new team members
* ───────────────────────────
* • V4.1 == "Adaptive" (OS-CFAR) mode. V2 == "Otsu" mode (py2 method_b).
* • Both detectors implement `WallDetector` so the same SweepInput drives both
* in parallel (drift-zero pairing).
* • BV is now from `BvEstimation.estimate()` (multi-method dispatch);
* sphere fit BV is kept as a cross-check value inside that result.
* • `gainDb` only affects V2; V4.1 ignores it on purpose.
* • All algorithms are pure functions — no side effects beyond logging.
* • For team docs see `docs/CHARLES-V41-API.md`,
* `docs/BV-CALCULATION-DESIGN.md`.
* ═════════════════════════════════════════════════════════════════════════
*/
package com.medithings.vesiscan.walldetect
import com.medithings.vesiscan.walldetect.algo.AnatomicalGate
import com.medithings.vesiscan.walldetect.algo.BModeScore
import com.medithings.vesiscan.walldetect.algo.BvEstimation
import com.medithings.vesiscan.walldetect.algo.BvFromSphere
import com.medithings.vesiscan.walldetect.algo.ChordConsensus
import com.medithings.vesiscan.walldetect.algo.Clipping
import com.medithings.vesiscan.walldetect.algo.ContrastAux
import com.medithings.vesiscan.walldetect.algo.DetectLumenFirst
import com.medithings.vesiscan.walldetect.algo.Geometry
import com.medithings.vesiscan.walldetect.algo.ImpulseReject
import com.medithings.vesiscan.walldetect.algo.PeakDetection
import com.medithings.vesiscan.walldetect.algo.SphereFit2Step
import com.medithings.vesiscan.walldetect.algo.SubsampleRefine
import com.medithings.vesiscan.walldetect.algo.WaveletDenoise
import com.medithings.vesiscan.walldetect.core.WdConfig
import com.medithings.vesiscan.walldetect.dto.ChannelResult
import com.medithings.vesiscan.walldetect.dto.DetectionResult
import com.medithings.vesiscan.walldetect.dto.DetectionSummary
import com.medithings.vesiscan.walldetect.dto.IntRange2
import com.medithings.vesiscan.walldetect.dto.ScoreSubscores
import com.medithings.vesiscan.walldetect.dto.SweepInput
import com.medithings.vesiscan.walldetect.dto.V41Diagnostics
import com.medithings.vesiscan.walldetect.dto.V41SphereFit
import com.medithings.vesiscan.walldetect.dto.Vec3
import com.medithings.vesiscan.walldetect.dto.WaveletCoefs
import com.medithings.vesiscan.walldetect.dto.WaveletEnergyRatio
import kotlin.math.sqrt
class V41Detector(
private val gatePreset: AnatomicalGate.Preset = AnatomicalGate.PHANTOM_530,
private val useLrTilt: Boolean = WdConfig.USE_LR_TILT
) : WallDetector {
override val algorithmId: String = DetectorIds.V4_1
// v4.1.1 (2026-04-29): three-regime validation upgrade.
// • Two-pass Hampel impulse rejection (ImpulseReject.detectWithLumenClean).
// • V41_PARAMS in DetectLumenFirst (mergeGapMax=5, gapPeakMargin=50,
// maxCandidatesPost=64) — phantom-on-rigid-floor speckle clusters
// and far-but-dominant wall+floor merged peaks now correctly handled.
// • V41_WEIGHTS in BModeScore (adds u_wamp + u_stl, Allen 1978 STA/LTA;
// phantom-mode trust no longer collapses with u_far→0).
// • PHANTOM_530 gate band [10, 60] mm — Neyman-Pearson loose prior.
// • New algo modules: StaLta, ChordConsensus (Tukey 1977 +
// Fischler-Bolles 1981), ImpulseReject, MorphClose (kept disabled).
// • Validated on 3-capture set: Center 1.5%, Corner 18%, 500 mL
// Phantom on Floor 0.3% BV error.
override val algorithmVersion: String = "v4.1.1"
override fun detect(input: SweepInput): DetectionResult {
val t0 = System.nanoTime()
val perCh: List<ChannelResult> = (0..5).map { ch -> detectChannel(ch, input) }
// ── 5/6. sweep-level sphere fit (Mode A only — Q5: gated < 4 → null) ──
val gatedDetections: List<Geometry.Detection> = perCh.map { cr ->
if (cr.v41Diag?.gated == true && cr.antIdx != null && cr.postIdx != null)
Geometry.Detection(cr.antIdx, cr.postIdx)
else
Geometry.Detection(null, null)
}
val gatedCount = gatedDetections.count { it.ant != null && it.post != null }
val v41Sphere: V41SphereFit? = if (gatedCount >= 4) {
val wallPts = Geometry.buildWallPointsVisual(gatedDetections, useLr = useLrTilt)
val fit = SphereFit2Step.fit2Step(
wallPts.map { it.xyz },
SphereFit2Step.Mode.AUTO
)
if (fit != null) buildSphereDto(fit, wallPts) else null
} else null
// ── BV dispatch (PR-13 — see docs/BV-CALCULATION-DESIGN.md) ──
// Use the gated ant/post pairs (post anatomical gate) as input. Sphere
// fit BV is passed as cross-check.
val bvDispatch = run {
// v4.1.1 — pass the B-mode score per channel so BvEstimation can
// run ChordConsensus (score-trust + Tukey/Fischler-Bolles).
// Score is only meaningful for gated detections; ungated channels
// get score=0 so the consensus filter excludes them anyway.
val dets = perCh.map { cr ->
val isGated = cr.v41Diag?.gated == true
BvEstimation.Detection(
ant = if (isGated) cr.antIdx else null,
post = if (isGated) cr.postIdx else null,
score = if (isGated) (cr.v41Diag?.score?.toDouble() ?: 0.0) else 0.0,
)
}
BvEstimation.estimate(
walls = dets,
sphereCrossCheckBvMl = v41Sphere?.bvMl,
)
}
// ── summary ──
val scores = perCh.mapNotNull { it.v41Diag?.score }
val chords = perCh.mapNotNull { it.chordMm }
val gated = perCh.count { it.v41Diag?.gated == true }
val chordMean: Float? = chords.takeIf { it.isNotEmpty() }?.average()?.toFloat()
val chordStd: Float? = stdF(chords)
val scoreMean: Float? = scores.takeIf { it.isNotEmpty() }?.let { it.average().toFloat() }
val scoreMin: Float? = scores.minOrNull()
val scoreMax: Float? = scores.maxOrNull()
val sweepTier: String? = scoreMean?.let { BModeScore.classify(it.toDouble()).tier }
return DetectionResult(
requestId = input.requestId,
timestampMs = input.timestampMs,
algorithm = algorithmId,
algorithmVersion = algorithmVersion,
processingMs = (System.nanoTime() - t0) / 1_000_000.0,
perChannel = perCh,
summary = DetectionSummary(
matchCount = perCh.count { it.antIdx != null && it.postIdx != null },
chordMmMean = chordMean,
chordMmStd = chordStd,
tier = sweepTier,
scoreMean = scoreMean,
scoreMin = scoreMin,
scoreMax = scoreMax,
gatedCount = gated,
v41Sphere = v41Sphere,
bvDispatch = bvDispatch,
)
)
}
private fun detectChannel(ch: Int, input: SweepInput): ChannelResult {
val rawIntList = input.adc[ch]
val rawInt = IntArray(rawIntList.size) { rawIntList[it] }
val raw = DoubleArray(rawInt.size) { rawInt[it].toDouble() }
// 1.b/c/d Wavelet (raw envelope)
val wlDenoised = WaveletDenoise.denoise(raw, levels = 3)
val wlDiag = WaveletDenoise.diagnose(raw, levels = 3)
val wlDecomp = WaveletDenoise.dwt(raw, levels = 3)
// 1.a + 2.a + 2.d + 3.a + 3.b — V4.1 two-pass detection.
// • Pass 1: DetectLumenFirst with V41_PARAMS (Hampel-aware
// mergeGapMax=5, gapPeakMargin=50, maxCandidatesPost=64).
// • Hampel-in-range: Hampel impulse rejection (Hampel 1974)
// restricted to the coarse lumen [ant+1, post-1] — handles
// 1-3 sample isolated impulses without touching wall samples.
// • Pass 2: re-detect on cleaned envelope. The two-pass
// structure is the architectural strength described in §6.5b.
val twoPass = ImpulseReject.detectWithLumenClean(
rawInt, DetectLumenFirst.Mode.Adaptive,
params = DetectLumenFirst.V41_PARAMS
)
val det = twoPass.refined
val sg = det.sg
val cfarThr = det.cfarThr ?: DoubleArray(sg.size) { det.adaptiveT } // safety
// 2.b clipping (raw)
val clipping = Clipping.isClipped(rawInt)
// 2.c peaks (sg, all candidates)
val peaksAll = PeakDetection.findPeaks1D(sg).toList()
// 3.c subsample refine (parabolic, peak kind) — operates on RAW envelope.
// SG smoothing slightly biases the parabola vertex, so refine uses the
// raw amplitude. The two-pass detector's `cleanedEnvelope` is used so
// that intra-lumen impulses (already removed in pass 2) don't pull the
// parabola during refinement.
val refineSubstrate = twoPass.cleanedEnvelope
val antRefined = det.ant?.let { SubsampleRefine.refineParabolic(refineSubstrate, it, SubsampleRefine.Kind.PEAK) }
val postRefined = det.post?.let { SubsampleRefine.refineParabolic(refineSubstrate, it, SubsampleRefine.Kind.PEAK) }
val antMmRaw: Float? = antRefined?.let { Geometry.sampleToMm(it).toFloat() }
?: det.ant?.let { Geometry.sampleToMm(it.toDouble()).toFloat() }
val postMmRaw: Float? = postRefined?.let { Geometry.sampleToMm(it).toFloat() }
?: det.post?.let { Geometry.sampleToMm(it.toDouble()).toFloat() }
// 4.a far-post contrast (sg)
val cR = ContrastAux.farPostContrast(sg, det.ant, det.post)
val sContrast: Float = (cR?.contrast ?: 0.0).toFloat()
val sContrastTier: String = ContrastAux.tierBand(cR?.contrast)
// 4.b B-mode composite score (RAW envelope) — V4.1 weighting.
// far 0.15 + dark 0.10 + ant 0.05 + post 0.05
// + wamp 0.45 (wall-peak amplitude vs lumen baseline)
// + stalta 0.20 (Allen-1978 STA/LTA impulse purity)
// The wamp + stalta pair handles the phantom-on-rigid-floor
// regime where there is no tissue echo behind the wall and the
// legacy s_contrast / u_far drops to 0.
val sR = BModeScore.score(
raw, det.ant, det.post,
weights = BModeScore.V41_WEIGHTS
)
val score: Float = (sR?.total ?: 0.0).toFloat()
val scoreSub: ScoreSubscores = if (sR != null) ScoreSubscores(
uFarPost = sR.sub.far.toFloat(),
uLumDark = sR.sub.dark.toFloat(),
uAntGrad = sR.sub.antGrad.toFloat(),
uPostGrad = sR.sub.postGrad.toFloat()
) else ScoreSubscores(0f, 0f, 0f, 0f)
val chTier: String = sR?.tier ?: "—"
// 4.c anatomical gate
val gate = AnatomicalGate.apply(
AnatomicalGate.Detection(det.ant, det.post),
channelIndex = ch,
preset = gatePreset
)
val gated = gate.passed
// 6.b single-channel BV (gated only)
val singleBv: Float? = if (gated && det.ant != null && det.post != null) {
BvFromSphere.bvFromChord((det.post - det.ant).toDouble()).toFloat()
} else null
// spans (DetectLumenFirst.detect already merged them)
val spansForDto: List<IntRange2> = det.spans.map { IntRange2(it.start, it.end) }
val chordMmFinal: Float? =
if (gated && antMmRaw != null && postMmRaw != null) postMmRaw - antMmRaw else null
// Wavelet diagnostic — diagnose() returns null for n < 8; sg.size = 100 → safe
val (rL1, rL2, rL3) = if (wlDiag != null) Triple(
wlDiag.ratio[0].toFloat(),
wlDiag.ratio[1].toFloat(),
wlDiag.ratio[2].toFloat()
) else Triple(0f, 0f, 0f)
return ChannelResult(
ch = ch,
sg = sg.map { it.toFloat() },
threshold = cfarThr.map { it.toFloat() },
antIdx = if (gated) det.ant else null,
postIdx = if (gated) det.post else null,
antRefined = if (gated) antRefined?.toFloat() else null,
postRefined = if (gated) postRefined?.toFloat() else null,
antMm = if (gated) antMmRaw else null,
postMm = if (gated) postMmRaw else null,
lumenStart = if (gated) det.lowStart else null,
lumenEnd = if (gated) det.lowEnd else null,
chordMm = chordMmFinal,
clipping = clipping,
v41Diag = V41Diagnostics(
waveletDenoised = wlDenoised.map { it.toFloat() },
waveletEnergyRatio = WaveletEnergyRatio(rL1, rL2, rL3),
waveletCoefs = WaveletCoefs(
l1 = wlDecomp.details[0].map { it.toFloat() },
l2 = wlDecomp.details[1].map { it.toFloat() },
l3 = wlDecomp.details[2].map { it.toFloat() },
a3 = wlDecomp.approx.map { it.toFloat() }
),
peaksAll = peaksAll,
spans = spansForDto,
score = score,
scoreSub = scoreSub,
gated = gated,
tier = chTier,
sContrast = sContrast,
sContrastTier = sContrastTier,
singleChannelBvMl = singleBv
)
)
}
/**
* Build V41SphereFit DTO from the SphereFit2Step result.
* Circle mode (constY/X/Z): inflate 2D centre back to 3D using the const-axis value.
* Sphere mode: 3D centre directly.
*/
private fun buildSphereDto(
fit: SphereFit2Step.FitResult,
wallPts: List<Geometry.WallPoint>
): V41SphereFit {
val center3 = if (fit.mode == "circle") {
val axes = fit.axes!!
val dropAxis = fit.dropAxis!!
val constVal = if (wallPts.isNotEmpty()) wallPts[0].xyz[dropAxis] else 0.0
val c = DoubleArray(3)
c[axes[0]] = fit.lmCenter[0]
c[axes[1]] = fit.lmCenter[1]
c[dropAxis] = constVal
c
} else {
fit.lmCenter
}
val rMm = fit.lmR
val bvMl = BvFromSphere.bvFromSphere(rMm)
return V41SphereFit(
mode = "A",
center = Vec3(center3[0].toFloat(), center3[1].toFloat(), center3[2].toFloat()),
radiusMm = rMm.toFloat(),
bvMl = bvMl.toFloat(),
residualStdMm = fit.residualStd.toFloat(),
wallPoints = wallPts.map {
Vec3(it.xyz[0].toFloat(), it.xyz[1].toFloat(), it.xyz[2].toFloat())
},
deltaRMm = (rMm - WdConfig.PHANTOM_530_R_MM).toFloat(),
deltaBvMl = (bvMl - WdConfig.PHANTOM_530_BV_ML).toFloat(),
nPoints = wallPts.size
)
}
private fun stdF(values: List<Float>): Float? {
if (values.isEmpty()) return null
if (values.size == 1) return 0f
val mu = values.average()
var sq = 0.0
for (v in values) { val d = v - mu; sq += d * d }
return sqrt(sq / values.size).toFloat()
}
}
@@ -0,0 +1,47 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* ─────────────────────────────────────────────────────────────────────────
* WallDetector — common contract for V2 and V4.1 detectors
* ─────────────────────────────────────────────────────────────────────────
*
* Two implementations exist:
* • V2Detector — py2 method_b (Otsu adaptive threshold + lumen-first)
* • V41Detector — Charles V4.1 (OS-CFAR + B-mode score + sphere fit)
*
* Both consume the same `SweepInput` and emit the same `DetectionResult`
* shape so a single sweep can drive both detectors in parallel
* (drift-zero pairing). Only the `algorithm` field differentiates output.
*
* Invariants (must hold for every implementation)
* ───────────
* • Single method: `detect(SweepInput): DetectionResult`
* • Idempotent: same input → same output (modulo `processingMs`)
* • No side effects beyond logging (no IO / no mutation of input)
* • Thread-safe: a single instance can be called from multiple coroutines
* • `gainDb`: V2 scales LOW_ECHO_AMP by 10^(-gainDb/20); V4.1 ignores it (Q4)
*/
package com.medithings.vesiscan.walldetect
import com.medithings.vesiscan.walldetect.dto.DetectionResult
import com.medithings.vesiscan.walldetect.dto.SweepInput
interface WallDetector {
/** Stable algorithm identifier. `"v2"` | `"v4_1"`. */
val algorithmId: String
/** Algorithm semantic version (bump when the numeric definition changes). */
val algorithmVersion: String
/**
* Run the detector on a single sweep. Returns a complete `DetectionResult`
* including per-channel diagnostics and a sweep-level summary.
*
* @param input 6 × 100 ADC matrix + probe metadata
* @return DetectionResult with `algorithm = algorithmId`
*/
fun detect(input: SweepInput): DetectionResult
}
@@ -0,0 +1,156 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/anatomical_gate.js (1:1).
*
* Biological-plausibility gate. Two priors:
* 1. Anterior-wall vertical depth band [antDepthMin, antDepthMax]
* probe-to-anterior depth, projected by cos(LR)·cos(SI), in band.
* PHANTOM_530 (BP2 530 mL): [22, 60] mm ← V41Detector default (Q6)
* CLINICAL (free bladder): [20, 70] mm
* 2. Wall-to-wall chord band (along beam) [chordMin, 2·R_max + chordSlack]
* PHANTOM_530: [5, 110.4] mm
* CLINICAL: [5, 130] mm
*
* Beam-angle correction uses Snell-refracted angles (PROBE.DEGREE / DEGREE_LR),
* NOT the mechanical CAD angles (MECH_DEGREE / MECH_DEGREE_LR).
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdConfig
import com.medithings.vesiscan.walldetect.core.WdProbe
import kotlin.math.cos
object AnatomicalGate {
private const val DEG = Math.PI / 180.0
data class Preset(
val antDepthMin: Double,
val antDepthMax: Double,
val rMax: Double,
val chordMin: Double,
val chordSlack: Double
)
/**
* Legacy in-vivo preset (kept for backwards compatibility).
* Tight ant-depth band — appropriate for human captures with the
* standard 12-25 mm abdominal-wall layer between probe and bladder.
*/
val LEGACY_INVIVO = Preset(
antDepthMin = 22.0, antDepthMax = 60.0,
rMax = 50.20, chordMin = 5.0, chordSlack = 10.0
)
/**
* V4.1 default — Neyman-Pearson "loose prior" for the 530 mL BP2
* phantom AND phantom-on-rigid-floor AND corner geometries. Tight
* discrimination (in-vivo vs reverberation) is delegated to the
* likelihood-ratio test (B-mode score, trust threshold 0.50). The
* gate only rejects what no acquisition geometry could ever produce.
* ant_vd > 10 mm — minimum probe near-field + coupling layer.
* ant_vd < 60 mm — extreme corner/off-axis still intersects bladder.
* chord ∈ [5, 110.4] — geometric chord of a sphere R ≤ 50.2 mm.
* Reference: Neyman J, Pearson ES. "On the problem of the most
* efficient tests of statistical hypotheses." Phil Trans R Soc A
* 231:289-337, 1933. Lehmann EL "Testing Statistical Hypotheses"
* 1986 §3 (loose prior + sharp likelihood for nuisance-parameter
* problems).
*/
val PHANTOM_530 = Preset(
antDepthMin = 10.0, antDepthMax = 60.0,
rMax = 50.20, chordMin = 5.0, chordSlack = 10.0
)
/** Free-bladder clinical preset. */
val CLINICAL = Preset(
antDepthMin = 20.0, antDepthMax = 70.0,
rMax = 65.0, chordMin = 5.0, chordSlack = 0.0
)
/** Per-channel detection input — only the indices matter for gating. */
data class Detection(val ant: Int?, val post: Int?)
data class BeforeGate(
val ant: Int,
val post: Int,
val antVdMm: Double,
val chordMm: Double
)
/**
* Gate decision. On pass: `ant`/`post` retained, `gateRejected = null`.
* On reject: `ant`/`post` set to null, `gateRejected` carries reason,
* `beforeGate` retains original indices + diagnostics.
*/
data class Result(
val ant: Int?,
val post: Int?,
val gateRejected: String?,
val antVdMm: Double?,
val chordMm: Double?,
val beforeGate: BeforeGate?
) {
val passed: Boolean get() = gateRejected == null && ant != null && post != null
}
fun apply(
detection: Detection,
channelIndex: Int,
preset: Preset = PHANTOM_530,
dps: Double = WdConfig.DPS_DEFAULT,
delay: Double = WdConfig.DELAY_MM_DEFAULT,
siAngles: DoubleArray = WdProbe.DEGREE,
lrAngles: DoubleArray = WdProbe.DEGREE_LR
): Result {
val ant = detection.ant
val post = detection.post
if (ant == null || post == null) {
return Result(
ant = null, post = null,
gateRejected = null,
antVdMm = null, chordMm = null, beforeGate = null
)
}
val si = siAngles.getOrElse(channelIndex) { 0.0 }
val lr = lrAngles.getOrElse(channelIndex) { 0.0 }
val cosBeamY = cos(lr * DEG) * cos(si * DEG)
val antPath = ant * dps + delay
val antVd = antPath * cosBeamY // vertical depth (mm)
val chord = (post - ant) * dps // along-beam (mm)
val chordMax = 2.0 * preset.rMax + preset.chordSlack
val rejected: String? = when {
antVd < preset.antDepthMin ->
"ant depth ${"%.1f".format(antVd)} mm < ${"%.0f".format(preset.antDepthMin)} mm"
antVd > preset.antDepthMax ->
"ant depth ${"%.1f".format(antVd)} mm > ${"%.0f".format(preset.antDepthMax)} mm"
chord < preset.chordMin ->
"chord ${"%.1f".format(chord)} mm < ${"%.0f".format(preset.chordMin)} mm"
chord > chordMax ->
"chord ${"%.1f".format(chord)} mm > ${"%.0f".format(chordMax)} mm " +
"(2·R_max + ${"%.0f".format(preset.chordSlack)} mm)"
else -> null
}
return if (rejected != null) {
Result(
ant = null, post = null,
gateRejected = rejected,
antVdMm = null, chordMm = null,
beforeGate = BeforeGate(ant, post, antVd, chord)
)
} else {
Result(
ant = ant, post = post,
gateRejected = null,
antVdMm = antVd, chordMm = chord, beforeGate = null
)
}
}
}
@@ -0,0 +1,230 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/bmode_score.js (1:1).
*
* Quantitative 0..1 B-mode confidence score per channel.
*
* Subscores (computed on RAW envelope, indices ant/post inclusive of walls,
* lumen STRICTLY between them):
*
* L = r[ant+1 .. post−1] (lumen)
* F = r[post+τ .. post+τ+W] (far-post window)
* O = r[0..ant] ∪ r[post..N] (outside lumen — JS strict: k > ant && k < post)
*
* s_far = max(0, mean(F) − mean(L)) far-post recovery (ADC)
* s_dark = max(0, mean(O) − mean(L)) lumen-darkness depth (ADC)
* g_ant = |r[ant+1] − r[ant−1]| / 2 anterior wall gradient (ADC)
* g_post = |r[post+1] − r[post−1]| / 2 posterior wall gradient (ADC)
*
* Soft saturation:
* u(x; x_50) = x / (x + x_50)
*
* Total:
* score = w_far · u_far + w_dark · u_dark + w_ant · u_ant + w_post · u_post
*
* Tier from score:
* ≥ 0.70 → high
* ≥ 0.40 → moderate
* ≥ 0.15 → low
* else → zero
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
import kotlin.math.min
object BModeScore {
// Defaults — DO NOT CHANGE without bumping algorithmVersion (golden tests will fail).
const val DEFAULT_TAU = 6
const val DEFAULT_WIN = 10
val DEFAULT_X50 = X50(far = 100.0, dark = 150.0, grad = 250.0, wamp = 200.0, stalta = 1.0)
val DEFAULT_WEIGHTS = Weights(far = 0.45, dark = 0.25, ant = 0.15, post = 0.15, wamp = 0.0, stalta = 0.0)
/**
* V4.1 preset (this work): adds u_wamp (wall-peak amplitude vs lumen
* baseline, Q6.5d) and u_stl (Allen-1978 STA/LTA impulse purity,
* StaLta.peakRatio − 1.0). Re-weights to make wamp the dominant term
* because:
* • u_far → 0 in phantom-on-rigid-floor regime where there is no
* tissue echo behind the wall;
* • integer ant/post often lands at the wall PEAK (gradient ≈ 0
* with neighbours), so u_ant/u_post are unreliable;
* • wall-peak amplitude is the most direct evidence of a real wall
* and survives intact across all geometry regimes.
* Weights sum to 1.0: far 0.15, dark 0.10, ant 0.05, post 0.05,
* wamp 0.45, stalta 0.20.
*/
val V41_WEIGHTS = Weights(
far = 0.15, dark = 0.10, ant = 0.05, post = 0.05, wamp = 0.45, stalta = 0.20
)
data class X50(
val far: Double,
val dark: Double,
val grad: Double,
val wamp: Double = 200.0,
val stalta: Double = 1.0
)
data class Weights(
val far: Double,
val dark: Double,
val ant: Double,
val post: Double,
val wamp: Double = 0.0,
val stalta: Double = 0.0
)
data class Subscores(
val far: Double,
val dark: Double,
val antGrad: Double,
val postGrad: Double,
val wallAmp: Double = 0.0,
val staLta: Double = 0.0
)
data class RawValues(
val sFar: Double,
val sDark: Double,
val gAnt: Double,
val gPost: Double,
val lumenMean: Double,
val farMean: Double?,
val outsideMean: Double,
val sWamp: Double = 0.0,
val rStl: Double = 0.0,
val pMax: Double = 0.0
)
data class Windows(
val lLo: Int, val lHi: Int, val fLo: Int, val fHi: Int,
val pLo: Int = 0, val pHi: Int = 0
)
data class Tier(val tier: String, val desc: String)
data class Result(
val total: Double,
val tier: String,
val desc: String,
val sub: Subscores,
val raw: RawValues,
val windows: Windows
)
/** Soft saturation u(x) = x / (x + x50); 0 if x or x50 ≤ 0. */
fun softSat(x: Double, x50: Double): Double {
if (x <= 0.0 || x50 <= 0.0) return 0.0
return x / (x + x50)
}
fun classify(score: Double?): Tier {
if (score == null) return Tier("—", "no detection")
return when {
score >= 0.70 -> Tier("high", "strong B-mode signature on all axes")
score >= 0.40 -> Tier("moderate", "B-mode signature attenuated but coherent")
score >= 0.15 -> Tier("low", "marginal — single-feature support")
else -> Tier("zero", "absent / shadowed / mis-detection")
}
}
/**
* Compute the composite B-mode score on RAW envelope.
* @param envelope raw envelope (NOT sg-smoothed; gradient subscores need sharp transitions)
* @return null if ant/post are invalid (null, equal, or pathological pair)
*/
fun score(
envelope: DoubleArray?,
ant: Int?,
post: Int?,
tau: Int = DEFAULT_TAU,
win: Int = DEFAULT_WIN,
x50: X50 = DEFAULT_X50,
weights: Weights = DEFAULT_WEIGHTS
): Result? {
if (envelope == null || envelope.isEmpty()) return null
if (ant == null || post == null) return null
if (ant >= post - 1) return null
val n = envelope.size
// Lumen mean (strictly between walls)
var lSum = 0.0
var lN = 0
for (k in (ant + 1) until post) { lSum += envelope[k]; lN++ }
val lumenMean = if (lN > 0) lSum / lN else 0.0
// Far-post window
val fLo = post + tau
val fHi = min(n, fLo + win)
var fSum = 0.0
var fN = 0
for (k in fLo until fHi) { fSum += envelope[k]; fN++ }
val farMean: Double? = if (fN > 0) fSum / fN else null
// Outside-lumen mean (everything except lumen interval)
// JS: `if (k > ant && k < post) continue;` — strict inside skipped, walls included
var oSum = 0.0
var oN = 0
for (k in 0 until n) {
if (k > ant && k < post) continue
oSum += envelope[k]; oN++
}
val outsideMean = if (oN > 0) oSum / oN else 0.0
// V4.1 — wall-peak window centred on detected post (8 samples).
// Used by u_wamp (wall-peak amplitude vs lumen baseline). A small
// bracket so a 1-sample mis-snap of post does not under-measure
// the wall height.
val pLo = maxOf(0, post - 2)
val pHi = minOf(n - 1, post + 5)
var pMax = Double.NEGATIVE_INFINITY
for (k in pLo..pHi) if (envelope[k] > pMax) pMax = envelope[k]
// V4.1 — STA/LTA impulse purity at the wall position (Allen 1978).
// Computed only when the wamp weight is non-zero (V4.1 preset);
// V2 / legacy presets skip this step entirely.
val rStl: Double = if (weights.stalta > 0.0)
maxOf(0.0, StaLta.peakRatio(envelope, post) - 1.0) else 0.0
// Raw subscore values (ADC)
val sFar = if (farMean != null) maxOf(0.0, farMean - lumenMean) else 0.0
val sDark = maxOf(0.0, outsideMean - lumenMean)
val gAnt = if (ant >= 1 && ant <= n - 2)
abs(envelope[ant + 1] - envelope[ant - 1]) / 2.0 else 0.0
val gPost = if (post >= 1 && post <= n - 2)
abs(envelope[post + 1] - envelope[post - 1]) / 2.0 else 0.0
val sWamp = maxOf(0.0, pMax - lumenMean)
// Mapped subscores ∈ [0,1]
val uFar = softSat(sFar, x50.far)
val uDark = softSat(sDark, x50.dark)
val uAnt = softSat(gAnt, x50.grad)
val uPost = softSat(gPost, x50.grad)
val uWamp = softSat(sWamp, x50.wamp)
val uStl = softSat(rStl, x50.stalta)
val total = weights.far * uFar +
weights.dark * uDark +
weights.ant * uAnt +
weights.post * uPost +
weights.wamp * uWamp +
weights.stalta * uStl
val cls = classify(total)
return Result(
total = total,
tier = cls.tier,
desc = cls.desc,
sub = Subscores(uFar, uDark, uAnt, uPost, uWamp, uStl),
raw = RawValues(sFar, sDark, gAnt, gPost, lumenMean, farMean, outsideMean,
sWamp = sWamp, rStl = rStl, pMax = pMax),
windows = Windows(lLo = ant + 1, lHi = post, fLo = fLo, fHi = fHi,
pLo = pLo, pHi = pHi)
)
}
}
@@ -0,0 +1,480 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Channel-combination BV dispatch — see docs/BV-CALCULATION-DESIGN.md.
*
* Implements 4-layer model (A/B/C/D) and 15-row dispatch table:
* A. FrustumNoLR / FrustumLR — Tanaka frustum + cap (py2 _bv_core 단순화)
* B. SphereLM — Kasa→LM (already in SphereFit2Step)
* C. Verathon — chord-as-diameter + lr_prior
* D. None — sentinel
*
* Simplifications vs py2 _bv_core:
* • Bottom & top caps both use hemispheres (hemisphere fallback) — no parabolic
* S(y) refinement. Acceptable accuracy ±5-10% for the live-compare use case.
* • LR ratio: single-chord only when one lateral; mean of both ratios when
* both laterals are present (simplified vs py2 3-chord offset-aware).
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdConfig
import com.medithings.vesiscan.walldetect.core.WdProbe
import com.medithings.vesiscan.walldetect.dto.BvDispatchResult
import kotlin.math.PI
import kotlin.math.abs
import kotlin.math.cbrt
import kotlin.math.cos
import kotlin.math.max
import kotlin.math.sin
import kotlin.math.sqrt
object BvEstimation {
private const val DEG = PI / 180.0
private const val LR_PRIOR = 1.2f
private const val LR_NO_DETECTION = 1.0f
private const val AREA_K = PI / 4.0 // (π/4) D²
/** Per-channel detection input — ant/post indices + the v4.1 B-mode score. */
data class Detection(val ant: Int?, val post: Int?, val score: Double = 0.0)
/**
* Main entry point. `walls` size must be 6. v4.1.1 dispatch:
*
* 1. ChordConsensus.filter (Tukey 1977 + Fischler-Bolles 1981) drops
* channels with B-mode score < 0.40 OR chord geometrically
* inconsistent with the multi-channel median. Returns the trusted
* set + median + MAD + leader.
* 2. Re-derive nC / nL / lrRatio from the TRUSTED set.
* 3. Dispatch:
* nC ≥ 2 → Frustum (length ∝ nC)
* nC = 1 + nL ≥ 1 → ChordMedian (consensus across center+lateral)
* nC = 1 → Verathon (chord-as-D × √lr_prior)
* nC = 0 + total ≥ 2 → ChordMedian (laterals only — diagnostic)
* else → None
* 4. Always carries trusted/rejected/leader/median/mad through the result.
*/
fun estimate(
walls: List<Detection>,
sphereCrossCheckBvMl: Float? = null,
dps: Double = WdConfig.DPS_DEFAULT,
delayMm: Double = WdConfig.DELAY_MM_DEFAULT,
siDeg: DoubleArray = WdProbe.DEGREE,
sensorZ: DoubleArray = WdProbe.SENSOR_Z,
trustScore: Double = ChordConsensus.DEFAULT_TRUST_SCORE,
chordTol: Double = ChordConsensus.DEFAULT_CHORD_TOL,
kMad: Double = ChordConsensus.DEFAULT_K_MAD,
): BvDispatchResult {
require(walls.size == 6) { "walls must have 6 entries (one per channel)" }
// ── 1) ChordConsensus filter (score-trust + Tukey/Fischler-Bolles) ──
val ccDetections = walls.map {
ChordConsensus.Detection(ant = it.ant, post = it.post, score = it.score)
}
val consensus = ChordConsensus.filter(
ccDetections, dps, trustScore = trustScore,
chordTol = chordTol, kMad = kMad,
)
val trusted = consensus.trusted
// ── 2) Re-derive center / lateral counts on TRUSTED set ──
val centerIdx = (0..3).filter {
it in trusted && walls[it].ant != null && walls[it].post != null
}
val lateralIdx = (4..5).filter {
it in trusted && walls[it].ant != null && walls[it].post != null
}
val nC = centerIdx.size
val nL = lateralIdx.size
// ── LR ratio uses ONLY trusted channels too ──
val lrRatio = computeLrRatio(walls, centerIdx, lateralIdx, dps, delayMm, siDeg)
// ── 3) Dispatch ──
val primary = when {
nC >= 2 -> frustum(walls, centerIdx, dps, delayMm,
siDeg, sensorZ, nC, nL, lrRatio,
sphereCrossCheckBvMl)
nC == 1 && nL >= 1 -> chordMedian(consensus, nC, nL, lrRatio,
sphereCrossCheckBvMl, "Mode B (1C + ${nL}L)")
nC == 1 -> verathon(walls[centerIdx[0]], centerIdx[0],
dps, siDeg, lrRatio, nC, nL, sphereCrossCheckBvMl)
(nC + nL) >= 2 -> chordMedian(consensus, nC, nL, lrRatio,
sphereCrossCheckBvMl, "lateral-only consensus")
else -> none(nC, nL, lrRatio,
if ((nC + nL) == 0) "no trusted detection (consensus filter)"
else "single lateral — no SI / chord-median basis")
}
// ── 4) Attach consensus diagnostics to the result ──
val withConsensus = primary.copy(
trustedChannels = consensus.trusted.sorted(),
rejectedChannels = consensus.rejected.map {
com.medithings.vesiscan.walldetect.dto.RejectedChannel(
it.ch, it.chord.toFloat(), it.reason,
)
},
consensusMedianChordMm = consensus.median?.toFloat(),
consensusMadMm = consensus.mad?.toFloat(),
leaderCh = consensus.leaderCh,
)
return applyAnatomicalBounds(withConsensus)
}
// ──────────────────────────────────────────────────────────
// Model B — ChordMedian (consensus-filtered chord-as-diameter)
//
// BV = (4/3)π·(median_chord/2)³. Used when frustum can't form (nC<2)
// but we still have ≥ 2 trusted detections that agree geometrically.
// ──────────────────────────────────────────────────────────
private fun chordMedian(
consensus: ChordConsensus.Result,
nC: Int, nL: Int,
lrRatio: Float,
sphereCrossCheck: Float?,
modeNote: String,
): BvDispatchResult {
val medChord = consensus.median
if (medChord == null || medChord <= 0.0) {
return none(nC, nL, lrRatio, "ChordMedian: empty consensus median")
}
val rMm = medChord / 2.0
val bvMl = (4.0 / 3.0) * PI * rMm * rMm * rMm / 1000.0
// Confidence: starts at 0.55 (above Verathon's 0.40, below Frustum's
// ≥ 0.70). +0.05 per trusted channel beyond 1, capped at 0.85.
val n = consensus.trusted.size
val confidence = (0.55f + 0.05f * (n - 1).coerceAtLeast(0)).coerceAtMost(0.85f)
val warnings = mutableListOf<String>("ChordMedian — $modeNote (n=$n)")
if (consensus.mad != null && consensus.mad > 5.0) {
warnings += "wide chord MAD (${"%.1f".format(consensus.mad)} mm) — geometry uncertain"
}
return BvDispatchResult(
bvMl = bvMl.toFloat(),
rMm = rMm.toFloat(),
method = "ChordMedian",
confidence = confidence,
nCenter = nC,
nLateral = nL,
lrRatio = lrRatio,
warnings = warnings,
sphereCrossCheckBvMl = sphereCrossCheck,
)
}
// ──────────────────────────────────────────────────────────
// Model A — Frustum (Tanaka)
// ──────────────────────────────────────────────────────────
private fun frustum(
walls: List<Detection>,
centerIdx: List<Int>,
dps: Double,
delayMm: Double,
siDeg: DoubleArray,
sensorZ: DoubleArray,
nC: Int,
nL: Int,
lrRatio: Float,
sphereCrossCheck: Float?,
): BvDispatchResult {
// 1. per-channel D, S, y
data class Cross(val ch: Int, val d_mm: Double, val S_mm2: Double, val y_mm: Double, val a: Double)
val xs = centerIdx.map { ch ->
val w = walls[ch]
val dAnt = delayMm + (w.ant!!).toDouble() * dps
val dPost = delayMm + (w.post!!).toDouble() * dps
val theta = siDeg[ch] * DEG
val L_raw = dPost - dAnt
val D = L_raw * cos(theta)
val S = AREA_K * D * D * lrRatio
val dMid = 0.5 * (dAnt + dPost)
val y = sensorZ[ch] + dMid * sin(theta)
Cross(ch, D, S, y, D * 0.5)
}.sortedBy { it.y_mm }
// 2. Frustum core (truncated cone integration between adjacent cross-sections)
var vCore = 0.0
for (i in 0 until xs.size - 1) {
val h = abs(xs[i + 1].y_mm - xs[i].y_mm)
val s1 = xs[i].S_mm2
val s2 = xs[i + 1].S_mm2
val v = (h / 3.0) * (s1 + s2 + sqrt(max(0.0, s1 * s2)))
vCore += v
}
// 3. Bottom & top caps —
// bottom: hemisphere (anterior dome of bladder)
// top : CONE (sigmoid/posterior tapers — ½ of hemisphere)
// Empirically tuned 2026-04-29: hemisphere on top was over-estimating
// BV by ~30-50% (e.g. centred 530mL phantom → 781mL). Cone is the
// simplified equivalent of py2 truncated-sphere top cap.
val rBot = xs.first().a
val rTop = xs.last().a
val vBot = (2.0 / 3.0) * PI * rBot * rBot * rBot // hemisphere
val vTop = (1.0 / 3.0) * PI * rTop * rTop * rTop // cone (½ hemi)
val bvMm3 = vCore + vBot + vTop
val bvMl = (bvMm3 / 1000.0).toFloat()
val method = if (nL > 0) "FrustumLR" else "FrustumNoLR"
val confidence = computeConfidence(method, nC, nL)
val rEq = cbrt(3.0 * bvMl * 1000.0 / (4.0 * PI)).toFloat()
return BvDispatchResult(
bvMl = bvMl,
rMm = rEq,
method = method,
confidence = confidence,
nCenter = nC,
nLateral = nL,
lrRatio = lrRatio,
sphereCrossCheckBvMl = sphereCrossCheck,
)
}
/** 2-channel cone fallback: frustum between the two + hemisphere on each end. */
private fun coneFallback(
walls: List<Detection>,
centerIdx: List<Int>,
dps: Double,
delayMm: Double,
siDeg: DoubleArray,
sensorZ: DoubleArray,
nC: Int,
nL: Int,
lrRatio: Float,
sphereCrossCheck: Float?,
): BvDispatchResult = frustum(walls, centerIdx, dps, delayMm, siDeg, sensorZ,
nC, nL, lrRatio, sphereCrossCheck)
.copy(method = "ConeFallback", confidence = 0.55f)
// ──────────────────────────────────────────────────────────
// Model B — Sphere fit on a point list (re-uses SphereFit2Step)
// ──────────────────────────────────────────────────────────
private fun sphereOnly(
walls: List<Detection>,
gatedIdx: List<Int>,
dps: Double,
delayMm: Double,
siDeg: DoubleArray,
sensorZ: DoubleArray,
nC: Int,
nL: Int,
lrRatio: Float,
sphereCrossCheck: Float?,
): BvDispatchResult {
// Build 3D wall points using visual coords (mirrors Geometry.wallIdxToXyzVisual
// but inlined here to avoid dependency on full WdProbe.SENSOR_X array).
val sx = WdProbe.SENSOR_X
val lr = WdProbe.DEGREE_LR
val pts = mutableListOf<DoubleArray>()
for (ch in gatedIdx) {
val w = walls[ch]
for (idx in listOf(w.ant!!, w.post!!)) {
val dist = delayMm + idx.toDouble() * dps
val tSi = siDeg[ch] * DEG
val tLr = lr[ch] * DEG
pts += doubleArrayOf(
sx[ch] + dist * sin(tLr) * cos(tSi),
dist * cos(tLr) * cos(tSi),
sensorZ[ch] + dist * sin(tSi),
)
}
}
val fit = SphereFit2Step.fit2Step(pts, SphereFit2Step.Mode.AUTO)
if (fit == null || fit.lmR <= 0) {
return none(nC, nL, lrRatio, "sphere fit failed")
}
val rMm = fit.lmR
val bvMl = ((4.0 / 3.0) * PI * rMm * rMm * rMm / 1000.0).toFloat()
val confidence = computeConfidence("SphereLM", nC, nL)
return BvDispatchResult(
bvMl = bvMl,
rMm = rMm.toFloat(),
method = "SphereLM",
confidence = confidence,
nCenter = nC,
nLateral = nL,
lrRatio = lrRatio,
sphereCrossCheckBvMl = sphereCrossCheck,
)
}
// ──────────────────────────────────────────────────────────
// Model C — Verathon-style single-channel chord
// ──────────────────────────────────────────────────────────
private fun verathon(
wall: Detection,
ch: Int,
dps: Double,
siDeg: DoubleArray,
lrRatio: Float,
nC: Int,
nL: Int,
sphereCrossCheck: Float?,
): BvDispatchResult {
val a = wall.ant!!.toDouble()
val p = wall.post!!.toDouble()
val theta = siDeg[ch] * DEG
val D = (p - a) * dps * cos(theta) // SI-corrected chord
val rSi = D / 2.0 // assumed great-circle radius
// Use measured lr_ratio if lateral was detected (auxiliary), else prior 1.2.
val lrFactor = if (lrRatio > 1.05f) lrRatio else LR_PRIOR
// Ellipsoid (a,b,c) ≈ (R_si, R_si·lr, R_si): BV = (4/3)π·a·b·c
// = sphere(R_si) × lrFactor — note: × lrFactor (not √lrFactor)
// Compared to py2 single-channel chord BV which uses simple sphere
// (lr_factor = 1) — we add lr scaling because the live use-case knows
// the bladder is LR>AP (Sun 2024).
val bvMl = ((4.0 / 3.0) * PI * rSi * rSi * rSi / 1000.0 * lrFactor).toFloat()
val warnings = mutableListOf<String>()
if (D < 25.0) warnings += "single CH chord too short (${"%.0f".format(D)}mm)"
if (D > 100.0) warnings += "single CH chord too long (${"%.0f".format(D)}mm)"
warnings += "1 channel only — recommend re-scan with more probe coverage"
return BvDispatchResult(
bvMl = bvMl,
rMm = rSi.toFloat(),
method = "Verathon",
confidence = if (D in 30.0..80.0) 0.35f else 0.20f,
nCenter = nC,
nLateral = nL,
lrRatio = lrFactor,
warnings = warnings,
sphereCrossCheckBvMl = sphereCrossCheck,
)
}
// ──────────────────────────────────────────────────────────
// Sentinel
// ──────────────────────────────────────────────────────────
private fun none(
nC: Int,
nL: Int,
lrRatio: Float,
reason: String,
): BvDispatchResult = BvDispatchResult(
bvMl = null, rMm = null, method = "None",
confidence = 0f, nCenter = nC, nLateral = nL, lrRatio = lrRatio,
warnings = listOf(reason),
)
// ──────────────────────────────────────────────────────────
// LR ratio — simplified single/two-chord (py2 compute_lr_ratio core)
// ──────────────────────────────────────────────────────────
private fun computeLrRatio(
walls: List<Detection>,
centerIdx: List<Int>,
lateralIdx: List<Int>,
dps: Double,
delayMm: Double,
siDeg: DoubleArray,
): Float {
if (lateralIdx.isEmpty() || centerIdx.isEmpty()) return LR_NO_DETECTION
// Center reference chord (D_center) — average of available center channels'
// SI-corrected chord lengths.
val dCenter = centerIdx.mapNotNull { ch ->
val w = walls[ch]
if (w.ant == null || w.post == null) null
else {
val theta = siDeg[ch] * DEG
val L = (w.post - w.ant) * dps
L * cos(theta)
}
}.takeIf { it.isNotEmpty() }?.average() ?: return LR_NO_DETECTION
if (dCenter <= 0) return 1.2f
// Per-lateral chord ratio
val ratios = lateralIdx.mapNotNull { ch ->
val w = walls[ch]
if (w.ant == null || w.post == null) return@mapNotNull null
val alpha = siDeg[ch] * DEG
val beta = WdProbe.DEGREE_LR[ch] * DEG
val P = cos(alpha) * cos(beta)
val L = (w.post - w.ant) * dps
val Dlat = L * P
val r = Dlat / dCenter
if (r <= 0 || r >= 1.0) null else r
}
if (ratios.isEmpty()) return 1.2f
val avgRatio = ratios.average()
// Single-chord b = |y_mid|/sqrt(1 - r²) is too dependent on probe placement;
// fall back to a simple inverse-ratio prior: lr_raw = 1 / r (clamped).
val lrRaw = (1.0 / avgRatio).coerceIn(1.0, 1.6).toFloat()
// Shrinkage toward 1.2 prior with confidence based on (1 - avg)
val confidence = ((1.0 - avgRatio) / 0.10).coerceIn(0.0, 1.0).toFloat()
return LR_PRIOR + (lrRaw - LR_PRIOR) * confidence
}
// ──────────────────────────────────────────────────────────
// Confidence model
// ──────────────────────────────────────────────────────────
private fun computeConfidence(method: String, nC: Int, nL: Int): Float = when (method) {
"FrustumLR" -> when {
nC >= 4 && nL == 2 -> 0.95f
nC >= 4 && nL == 1 -> 0.85f
nC == 3 && nL == 2 -> 0.85f
nC == 3 && nL == 1 -> 0.75f
nC == 2 && nL >= 1 -> 0.55f
else -> 0.50f
}
"FrustumNoLR" -> when (nC) {
4 -> 0.80f; 3 -> 0.70f; 2 -> 0.55f; else -> 0.50f
}
"SphereLM" -> when {
nC + nL >= 4 -> 0.65f
nC == 1 && nL == 2 -> 0.55f
nC == 1 && nL == 1 -> 0.45f
nC == 0 && nL == 2 -> 0.30f
else -> 0.40f
}
"ConeFallback" -> 0.55f
"Verathon" -> 0.35f
else -> 0f
}
// ──────────────────────────────────────────────────────────
// Anatomical safety rails + cross-check
// ──────────────────────────────────────────────────────────
private fun applyAnatomicalBounds(r: BvDispatchResult): BvDispatchResult {
val bv = r.bvMl ?: return r
val warnings = r.warnings.toMutableList()
var conf = r.confidence
when {
bv < 5f -> warnings += "BV < 5 mL — likely empty / not detected"
bv < 30f -> warnings += "BV < 30 mL — low / verify"
bv > 600f -> warnings += "BV > 600 mL — verify probe placement"
bv > 1000f -> {
return r.copy(bvMl = null, method = "None", confidence = 0f,
warnings = warnings + "BV > 1000 mL — rejected (unrealistic)")
}
}
// Cross-check vs sphere if both available
r.sphereCrossCheckBvMl?.let { spBv ->
if (r.method.startsWith("Frustum")) {
val diff = abs(bv - spBv) / bv
if (diff > 0.15f) {
warnings += "method disagreement (frustum=%.0f vs sphere=%.0f)".format(bv, spBv)
conf *= 0.7f
}
}
}
return r.copy(warnings = warnings, confidence = conf)
}
}
@@ -0,0 +1,28 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/bv_from_sphere.js (1:1).
* V4 canonical BV from sphere radius (mm) → mL.
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdConfig
object BvFromSphere {
/** BV (mL) from sphere radius (mm). */
fun bvFromSphere(rMm: Double): Double =
(4.0 / 3.0) * Math.PI * rMm * rMm * rMm / 1000.0
/**
* Legacy: chord-as-diameter sphere (single-channel BV).
* @param chordSamples post − ant in samples
* @param dpsMm distance per sample (mm), default WdConfig.DPS_DEFAULT
*/
fun bvFromChord(chordSamples: Double, dpsMm: Double = WdConfig.DPS_DEFAULT): Double {
val chordMm = chordSamples * dpsMm
return bvFromSphere(chordMm / 2.0)
}
}
@@ -0,0 +1,183 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/chord_consensus.js (1:1).
*
* Multi-channel chord-consensus outlier rejection — V4.1 ONLY.
*
* References
* Tukey JW. "Exploratory Data Analysis." Addison-Wesley, 1977.
* Boxplot / IQR definition of "isolated outliers" (k·MAD threshold).
* Fischler MA, Bolles RC. "Random sample consensus." Comm ACM
* 24(6):381-395, 1981. doi:10.1145/358669.358692
* Leader-driven consensus paradigm used here for the N=2 case.
* Rousseeuw PJ, Croux C. "Alternatives to the median absolute
* deviation." J Am Stat Assoc 88(424):1273-1283, 1993.
* MAD with 1.4826 normalisation for asymptotic Gaussian consistency.
*
* Algorithm
* 1. Collect score-passing channels (score ≥ TRUST_SCORE).
* 2. Compute chord_mm = (post − ant) · dps for each.
* 3. CASE A (N=0): nothing trusted.
* CASE B (N=1): trust the only channel.
* CASE C (N=2): pair test — keep both unless
* |chord_a − chord_b| / max > τ_chord (default 0.25); on inconsistency
* keep only the higher-score channel (Fischler-Bolles leader-driven).
* CASE D (N≥3): Tukey isolated-outlier rejection — drop channels with
* |chord_i − median| > k_mad · 1.4826 · MAD (default k_mad = 2.0).
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.PI
import kotlin.math.abs
import kotlin.math.cbrt
import kotlin.math.max
object ChordConsensus {
const val DEFAULT_TRUST_SCORE = 0.40
/**
* Pair-test tolerance — chord deviation between two trusted channels.
*
* Sphere-geometry rationale: the Fischler-Bolles 1981 RANSAC pair
* test assumes HOMOGENEOUS measurements (multiple noisy estimates
* of the same value). Ultrasound bladder chords are NOT homogeneous:
* each beam crosses the sphere at a different offset d from the
* centre, yielding chord c = 2·√(R² − d²) varying naturally in
* [0, 2R]. For two "useful" chords (both ≥ R, half-diameter
* coverage) Euclidean geometry permits up to 50% pair deviation.
* The earlier 25% default rejected legitimate off-axis observations
* (Corner 530 CH2/CH3: 36.8% deviation, but both chords consistent
* with R = 50.20 mm sphere at d = 40 mm and d = 16 mm).
*/
const val DEFAULT_CHORD_TOL = 0.50
/**
* Tukey kMad scale on the N ≥ 3 MAD-based outlier threshold.
* Tukey's classical 2σ recommendation applies to homogeneous
* samples; for sphere chords we accept up to 4σ on the natural
* beam-offset distribution (same sphere-geometry rationale as
* DEFAULT_CHORD_TOL).
*/
const val DEFAULT_K_MAD = 4.0
data class Detection(
val ant: Int?,
val post: Int?,
val score: Double = 0.0
)
data class Rejected(val ch: Int, val chord: Double, val reason: String)
data class Result(
val trusted: Set<Int>,
val median: Double?,
val mad: Double?,
val leaderCh: Int?,
val rejected: List<Rejected>
)
private fun median(arr: List<Double>): Double? {
if (arr.isEmpty()) return null
val s = arr.sorted()
return s[s.size / 2]
}
private fun mad(arr: List<Double>, med: Double): Double {
if (arr.isEmpty()) return 0.0
val dev = arr.map { abs(it - med) }
return median(dev) ?: 0.0
}
private data class Passer(val ch: Int, val score: Double, val chord: Double)
/**
* Run the consensus filter. Detections array is per-channel (size 6 typical).
* Channels with ant==null OR post==null OR score<trustScore are excluded.
*/
fun filter(
detections: List<Detection>,
dps: Double,
trustScore: Double = DEFAULT_TRUST_SCORE,
chordTol: Double = DEFAULT_CHORD_TOL,
kMad: Double = DEFAULT_K_MAD
): Result {
val passers = mutableListOf<Passer>()
detections.forEachIndexed { ch, r ->
if (r.ant == null || r.post == null) return@forEachIndexed
if (r.score < trustScore) return@forEachIndexed
passers.add(Passer(ch, r.score, (r.post - r.ant) * dps))
}
val trusted = mutableSetOf<Int>()
val rejected = mutableListOf<Rejected>()
var med: Double? = null
var madVal: Double? = null
var leader: Int? = null
if (passers.isEmpty()) return Result(trusted, null, null, null, rejected)
if (passers.size == 1) {
val p = passers[0]
trusted.add(p.ch)
return Result(trusted, p.chord, null, p.ch, rejected)
}
if (passers.size == 2) {
val sorted = passers.sortedByDescending { it.score }
val a = sorted[0]
val b = sorted[1]
val dev = abs(a.chord - b.chord) / max(a.chord, b.chord)
leader = a.ch
if (dev > chordTol) {
trusted.add(a.ch)
rejected.add(Rejected(b.ch, b.chord,
"pair-inconsistent (Δ ${"%.0f".format(dev * 100)}% > ${"%.0f".format(chordTol * 100)}%)"))
} else {
trusted.add(a.ch); trusted.add(b.ch)
}
med = median(passers.map { it.chord })
return Result(trusted, med, null, leader, rejected)
}
// N ≥ 3 — Tukey isolated-outlier rejection
val chords = passers.map { it.chord }
med = median(chords)!!
madVal = mad(chords, med)
val sigma = 1.4826 * madVal
val threshold = kMad * sigma
leader = passers.maxByOrNull { it.score }!!.ch
for (p in passers) {
val d = abs(p.chord - med)
if (madVal == 0.0 || d <= threshold) {
trusted.add(p.ch)
} else {
rejected.add(Rejected(p.ch, p.chord,
"isolated-outlier (|Δ| ${"%.1f".format(d)} mm > ${"%.1f".format(threshold)} mm = ${kMad}·MAD)"))
}
}
return Result(trusted, med, madVal, leader, rejected)
}
data class MedianBV(val bvMl: Double?, val rMm: Double?, val n: Int, val medianChord: Double?)
/**
* Tukey-robust chord-median BV on the trusted set: BV = (4/3)π(c/2)³.
*/
fun medianBV(detections: List<Detection>, dps: Double, trusted: Set<Int>): MedianBV {
val chords = mutableListOf<Double>()
detections.forEachIndexed { ch, r ->
if (ch in trusted && r.ant != null && r.post != null) {
chords.add((r.post - r.ant) * dps)
}
}
if (chords.isEmpty()) return MedianBV(null, null, 0, null)
val med = median(chords)!!
val rMm = med / 2.0
val bv = (4.0 / 3.0) * PI * rMm * rMm * rMm / 1000.0
return MedianBV(bv, rMm, chords.size, med)
}
}
@@ -0,0 +1,61 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/clipping.js (1:1).
* Dead-zone projection ψ_T(z)=max(z,T) + log envelope (dB).
* Reference: Donoho (1995) soft-thresholding; py4/clipping.py.
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdNumeric
import kotlin.math.log10
import kotlin.math.max
object Clipping {
/** Dead-zone projection ψ_T(z) = max(z, T) for scalar threshold. */
fun signalClip(sg: DoubleArray, threshold: Double): DoubleArray {
val n = sg.size
val out = DoubleArray(n)
for (i in 0 until n) out[i] = max(sg[i], threshold)
return out
}
/** Dead-zone projection with per-sample threshold. */
fun signalClip(sg: DoubleArray, threshold: DoubleArray): DoubleArray {
val n = sg.size
val out = DoubleArray(n)
for (i in 0 until n) out[i] = max(sg[i], threshold[i])
return out
}
/**
* Log envelope in dB: 20·log10(max(sg, eps) / base).
* If base is null, use max(median(sg), 1.0).
*/
fun logEnvelope(sg: DoubleArray, base: Double? = null, eps: Double = 1.0): DoubleArray {
val n = sg.size
val b = base ?: max(WdNumeric.median(sg), 1.0)
val out = DoubleArray(n)
for (i in 0 until n) out[i] = 20.0 * log10(max(sg[i], eps) / b)
return out
}
/**
* Detect ADC saturation: returns true if any sample exceeds the saturation
* threshold (default 4090 for 12-bit ADC).
* Used by V4.1 to flag channels where the wall echo is clipping the rail.
*/
fun isClipped(raw: DoubleArray, saturationThreshold: Double = 4090.0): Boolean {
for (v in raw) if (v >= saturationThreshold) return true
return false
}
/** Same, but on IntArray (typical sweep input). */
fun isClipped(raw: IntArray, saturationThreshold: Int = 4090): Boolean {
for (v in raw) if (v >= saturationThreshold) return true
return false
}
}
@@ -0,0 +1,112 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/contrast_aux.js (1:1).
*
* B-mode far-post contrast (CharlesKWON V4 reference, retained as a confidence
* metric in V4.1):
*
* Score = max(0, mean(sg[bw+τ : bw+τ+W]) − mean(sg[ant+1 : bw−1]))
*
* Tier (legacy, drives MultiModeBv trust criterion at τ=80):
* ≥ 200 → high
* ≥ 80 → moderate
* > 0 → low
* else → zero
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdConfig
import kotlin.math.max
import kotlin.math.min
object ContrastAux {
const val DEFAULT_TAU = 6
const val DEFAULT_WIN = 10
data class Result(
val contrast: Double,
val lumenMean: Double?,
val farMean: Double?,
val farLo: Int,
val farHi: Int,
val status: String // "positive" | "zero" | "truncated"
)
data class Tier(
val tier: String, // "high" | "moderate" | "low" | "zero" | "—"
val desc: String
)
fun farPostContrast(
sg: DoubleArray,
fw: Int?,
bw: Int?,
tau: Int = WdConfig.CONTRAST_TAU,
win: Int = WdConfig.CONTRAST_WIN
): Result? {
if (fw == null || bw == null) return null
val n = sg.size
val farLo = bw + tau
val farHi = min(n, farLo + win)
if (farLo >= farHi) {
return Result(
contrast = 0.0,
lumenMean = null,
farMean = null,
farLo = farLo,
farHi = farHi,
status = "truncated"
)
}
var farSum = 0.0
for (k in farLo until farHi) farSum += sg[k]
val farMean = farSum / (farHi - farLo)
val lumenMean: Double = if (bw - fw <= 1) {
sg[fw]
} else {
var s = 0.0
for (k in (fw + 1) until bw) s += sg[k]
s / (bw - fw - 1)
}
val contrast = max(0.0, farMean - lumenMean)
return Result(
contrast = contrast,
lumenMean = lumenMean,
farMean = farMean,
farLo = farLo,
farHi = farHi,
status = if (contrast > 0.0) "positive" else "zero"
)
}
/**
* Confidence label from contrast magnitude (legacy ADC scale).
* @param contrast null → "—" (no detection)
*/
fun classify(contrast: Double?): Tier {
if (contrast == null) return Tier("—", "no detection")
return when {
contrast >= 200 -> Tier("high", "canonical lumen-dark / far-post-bright pattern")
contrast >= 80 -> Tier("moderate", "far-post recovery present but attenuated")
contrast > 0 -> Tier("low", "marginal far-post brightness — verify")
else -> Tier("zero", "no brightness recovery — possible shadow / mis-detection")
}
}
/** Convenience tier string for V41Diagnostics.sContrastTier ("below80"|"80-200"|">=200"). */
fun tierBand(contrast: Double?): String {
if (contrast == null) return "below80"
return when {
contrast >= 200 -> ">=200"
contrast >= 80 -> "80-200"
else -> "below80"
}
}
}
@@ -0,0 +1,50 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/denoising.js (1:1).
* Savitzky-Golay (5,2) FIR with polynomial edge interpolation.
* Mirror of py4/denoising.sg_smooth (and py2/denoising.sg_smooth — 1:1 numerically).
*/
package com.medithings.vesiscan.walldetect.algo
object Denoising {
private val INNER = doubleArrayOf(-3.0, 12.0, 17.0, 12.0, -3.0).map { it / 35.0 }.toDoubleArray()
private val EDGE_LEFT_0 = doubleArrayOf(31.0, 9.0, -3.0, -5.0, 3.0).map { it / 35.0 }.toDoubleArray()
private val EDGE_LEFT_1 = doubleArrayOf(9.0, 13.0, 12.0, 6.0, -5.0).map { it / 35.0 }.toDoubleArray()
private val EDGE_RIGHT_1 = doubleArrayOf(-5.0, 6.0, 12.0, 13.0, 9.0).map { it / 35.0 }.toDoubleArray()
private val EDGE_RIGHT_0 = doubleArrayOf(3.0, -5.0, -3.0, 9.0, 31.0).map { it / 35.0 }.toDoubleArray()
/**
* Savitzky–Golay (window=5, poly=2) smoothing.
* Inner samples use the standard 5-tap kernel; edge samples (0/1/n-2/n-1)
* use polynomial-interpolation kernels matching numpy.polynomial fit.
*/
fun sgSmooth(x: DoubleArray): DoubleArray {
val n = x.size
if (n < 5) return x.copyOf()
val out = DoubleArray(n)
for (i in 2 until n - 2) {
var s = 0.0
for (k in 0 until 5) s += INNER[k] * x[i - 2 + k]
out[i] = s
}
var s0 = 0.0
var s1 = 0.0
var sn2 = 0.0
var sn1 = 0.0
for (k in 0 until 5) {
s0 += EDGE_LEFT_0[k] * x[k]
s1 += EDGE_LEFT_1[k] * x[k]
sn2 += EDGE_RIGHT_1[k] * x[n - 5 + k]
sn1 += EDGE_RIGHT_0[k] * x[n - 5 + k]
}
out[0] = s0
out[1] = s1
out[n - 2] = sn2
out[n - 1] = sn1
return out
}
}
@@ -0,0 +1,364 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/detect_lumen_first.js (1:1).
*
* V2 main detector (lumen-first wall pair).
* modes:
* Mode.Fixed(thr) → V2: lumen mask = sg ≤ thr (default 1150)
* Mode.Adaptive → V4.1: lumen mask = sg ≤ median(OS-CFAR(sg)) [PR-4]
* Mode.Scalar(value) → arbitrary scalar threshold (used by §3 7-method comparison)
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdNumeric
object DetectLumenFirst {
/**
* Default detection parameters — DO NOT CHANGE without bumping
* algorithmVersion. V2 uses these. V4.1 passes a different
* `ParamSet` via the new detect(raw, mode, params) overload below.
*/
object Params {
// Aligned to py2/for_app_share/config_6ch.py — LOW_ECHO_AMP bumped 1150 → 1250.
const val LOW_ECHO_AMP = 1250.0
const val LOW_MIN_LEN = 3
const val MERGE_GAP_MAX = 3
const val PEAK_SEARCH_WIN = 20
const val POST_MAX_IDX = 80
// Anatomical near-field cutoff (V4.1 ANT side). Sample 7 ≈
// half of the strict Fresnel near-field (Macovski 1983 §4:
// N = D²/4λ ≈ 35 mm ≈ sample 18 for the TB370FU 6 MHz / 6 mm
// aperture probe). Provides a structural margin against the
// steepest-angle channel (CH3, cos 0.9385) where the geometric
// near-field reaches deepest. See §6.5f of the verify paper.
const val ANT_MIN_IDX = 7
const val MIN_PEAK_MARGIN = 30.0
const val MIN_URINE_LEN = 3
const val GAP_PEAK_MARGIN = 5.0
// §3.7 (added 2026-05) — Stage 2 bimodal-merge dual-threshold
// guard ceiling, expressed as ADC margin above CFAR T. Wall+
// container double-peak failure mode (Japan-standard 150 mL
// CH0/CH5): a strong reflector beyond the bladder pulls Otsu's
// threshold up, mis-classifying the legitimate intermediate
// wall echo as speckle. The CFAR-derived ceiling (T + 250) is
// calibrated to lumen-noise statistics; the AND combination
// with otsuThr provides cross-validation across two orthogonal
// histograms. See SpanUtils.mergeBimodal §3.7 docstring.
const val GAP_PEAK_MARGIN_STAGE2 = 250.0
const val EDGE_DIST_DECAY = 0.12
const val VALLEY_STOP_RISE = 50.0
const val MAX_PEAK_CANDIDATES_ANT = 3
const val MAX_PEAK_CANDIDATES_POST = 3
}
/**
* Per-call parameter override. V4.1 uses V41_PARAMS; V2 leaves the
* argument null and inherits the static defaults.
*
* V4.1-specific changes from defaults:
* • MERGE_GAP_MAX 3 → 5
* lumen-first now closes gaps up to 5 contiguous "above-T"
* samples (phantom-on-rigid-floor speckle clusters are 4-5
* samples wide).
* • GAP_PEAK_MARGIN 5 → 50
* gap-peak guard requires a peak ≥ T+50 ADC to refuse the
* merge. Empirical separation: phantom speckle peaks 30-70
* ADC above T, real wall echoes 200-600 ADC above T.
* • MAX_PEAK_CANDIDATES_POST 3 → 64
* prominence comparison considers all candidates in the
* search window. Recovers the wall+floor merged echo on
* phantom captures where its distance rank exceeds 3.
*
* References for the parameter rationale:
* • Mathematical morphology — Serra 1982; Soille 2003.
* • Hampel impulse rejection (companion filter for narrow
* impulses ≤ 3 samples) — Hampel 1974; Pearson 2016.
*/
data class ParamSet(
val lowEchoAmp: Double = Params.LOW_ECHO_AMP,
val lowMinLen: Int = Params.LOW_MIN_LEN,
val mergeGapMax: Int = Params.MERGE_GAP_MAX,
val peakSearchWin: Int = Params.PEAK_SEARCH_WIN,
val postMaxIdx: Int = Params.POST_MAX_IDX,
val minPeakMargin: Double = Params.MIN_PEAK_MARGIN,
val minUrineLen: Int = Params.MIN_URINE_LEN,
val gapPeakMargin: Double = Params.GAP_PEAK_MARGIN,
/**
* §3.7 Stage 2 dual-threshold guard ceiling (V4.1 only). When
* positive AND `mode == Adaptive`, Stage 2 bimodal merge
* additionally requires the gap-peak amplitude to fall below
* `T + gapPeakMarginStage2`. AND-combined with the Otsu split
* to reject the wall+container double-peak pattern. 0 disables.
*/
val gapPeakMarginStage2: Double = 0.0,
val edgeDistDecay: Double = Params.EDGE_DIST_DECAY,
val valleyStopRise: Double = Params.VALLEY_STOP_RISE,
val maxPeakCandidatesAnt: Int = Params.MAX_PEAK_CANDIDATES_ANT,
val maxPeakCandidatesPost: Int = Params.MAX_PEAK_CANDIDATES_POST,
/**
* V4.1 anatomical near-field ANT floor (sample index). 0 disables
* the edge-fallback. See [Params.ANT_MIN_IDX] for the Macovski
* 1983 / Kremkau 2017 grounding.
*/
val antMinIdx: Int = 0,
/**
* Running-median pre-filter window (V4.1 only).
* 0 = disabled (V2 path stays bit-identical). 5 = V4.1 default,
* absorbs 1–2-sample isolated speckle bumps inside the lumen
* before OS-CFAR thresholding. See [MedianFilter] header for the
* Tukey 1974 / Justusson 1981 / Davies-Gather 1993 rationale.
*/
val medianWin: Int = 0,
)
/**
* V4.1 default parameter set — this work.
* Calibrated against the 3-capture validation set (Center 530 mL,
* Corner 530 mL, 500 mL Phantom on Floor) so all three regimes pass.
*/
val V41_PARAMS = ParamSet(
mergeGapMax = 5,
gapPeakMargin = 50.0,
gapPeakMarginStage2 = Params.GAP_PEAK_MARGIN_STAGE2,
maxPeakCandidatesPost = WallSelect.MAX_PEAK_CANDIDATES_POST,
medianWin = 7,
antMinIdx = Params.ANT_MIN_IDX,
)
/** Static defaults wrapped as a ParamSet. V2 path. */
val DEFAULT_PARAMS = ParamSet()
sealed interface Mode {
/** Fixed amplitude threshold (legacy V2). */
data class Fixed(val thr: Double = Params.LOW_ECHO_AMP) : Mode
/** V4.1 OS-CFAR adaptive (median of per-sample threshold). */
data object Adaptive : Mode
/** py2 method_b — 1-D Otsu over sg[0..POST_MAX_IDX]. */
data object Otsu : Mode
data class Scalar(val value: Double) : Mode
}
/**
* Detection result. Mirrors JS object literal returned by detect_lumen_first.detect().
* `cfarThr` is null in fixed/scalar modes.
*/
data class Result(
val mode: Mode,
val raw: DoubleArray,
val sg: DoubleArray,
val cfarThr: DoubleArray?,
val adaptiveT: Double,
val lowMask: BooleanArray,
val rawSpans: List<SpanUtils.Span>,
val spans: List<SpanUtils.Span>,
val ant: Int?,
val post: Int?,
val lowStart: Int?,
val lowEnd: Int?,
val lowMean: Double?,
val peakMin: Double?,
val urineLen: Int?,
val failReason: String?
)
fun detect(rawAdc: IntArray, mode: Mode = Mode.Fixed()): Result {
val raw = DoubleArray(rawAdc.size) { rawAdc[it].toDouble() }
return detect(raw, mode, DEFAULT_PARAMS)
}
fun detect(rawAdc: IntArray, mode: Mode, params: ParamSet): Result {
val raw = DoubleArray(rawAdc.size) { rawAdc[it].toDouble() }
return detect(raw, mode, params)
}
fun detect(rawAdc: DoubleArray, mode: Mode = Mode.Fixed()): Result =
detect(rawAdc, mode, DEFAULT_PARAMS)
fun detect(rawAdc: DoubleArray, mode: Mode, params: ParamSet): Result {
val raw = rawAdc.copyOf()
val sg = Denoising.sgSmooth(raw)
val n = sg.size
// V4.1: running median (Tukey 1974 / Justusson 1981) over the
// SG envelope BEFORE OS-CFAR + lumen-mask. Applied only to the
// adaptive path; the original `sg` is retained for wall-prominence
// (peak sharpness preserved). V2 / Otsu / Scalar leave
// sgForMask === sg → bit-identical to pre-filter behaviour.
val sgForMask: DoubleArray = if (mode is Mode.Adaptive && params.medianWin > 1) {
// Iterated to fixed-point (Justusson 1981 §4). Single pass
// leaves residual bumps in dense alternating clusters; 2
// iterations converge to the root signal.
MedianFilter.runningMedianRoot(sg, params.medianWin)
} else sg
val cfarThr: DoubleArray?
val T: Double
when (mode) {
is Mode.Adaptive -> {
// V4.1 OS-CFAR (Rohling 1983).
cfarThr = ThresholdOsCfar.perSample(sgForMask)
T = WdNumeric.median(cfarThr)
}
is Mode.Otsu -> {
// py2 method_b — Otsu over sg[0..POST_MAX_IDX] (skip far-tail).
cfarThr = null
val cap = minOf(params.postMaxIdx + 1, sg.size)
val view = DoubleArray(cap) { sg[it] }
T = Otsu.otsu1d(view)
}
is Mode.Scalar -> {
cfarThr = null
T = mode.value
}
is Mode.Fixed -> {
cfarThr = null
T = mode.thr
}
}
val lowMask = BooleanArray(n) { sgForMask[it] <= T }
val rawSpans = SpanUtils.contiguousTrueSpans(lowMask)
.filter { (it.end - it.start + 1) >= params.lowMinLen }
// Stage 1 — width-bounded amplitude-aware merge.
// Gap-peak amplitude is read from the ORIGINAL `sg` (not
// sgForMask) because the running median can clip a real wall
// peak to its plateau-median value, which on borderline cases
// drops below T + gapPeakMargin and would erroneously merge
// across the wall. The unfiltered sg preserves the true peak.
val gapPeakThr = T + params.gapPeakMargin
var spans = SpanUtils.mergeCloseSpans(rawSpans, params.mergeGapMax, sg, gapPeakThr)
// Stage 2 — V4.1 adaptive only. Otsu 1979 bimodal split + the
// anatomical postMaxIdx ceiling closes the wide-cluster case
// (e.g. 6-sample alternating speckle on the 150 mL Japan body
// phantom CH1) that exceeds the running median's ⌊W/2⌋ = 3
// absorption width.
//
// The gate `spans.size >= 3` restricts stage 2 to lumens
// that were FRAGMENTED by stage 1 — a normal capture leaves
// stage 1 with exactly 2 spans (lumen + post-wall tail) and
// needs no further merging. ≥ 3 spans signals an intra-lumen
// speckle cluster broke the lumen into pieces; only then is
// bimodal merging applied.
//
// The 1-ADC-resolution otsu1dInteger is used because the
// coarse 64-bin variant shifts the bimodal boundary by 10–20
// ADC and can flip the decision on borderline walls.
if (mode is Mode.Adaptive && spans.size >= 3) {
val cap = minOf(params.postMaxIdx + 1, sg.size)
val view = DoubleArray(cap) { sg[it] }
val otsuThr = Otsu.otsu1dInteger(view)
// §3.7 — pass T + GAP_PEAK_MARGIN_STAGE2 as the CFAR-derived
// dual-guard ceiling (null when disabled, preserving legacy
// single-threshold behaviour for non-V4.1 callers).
val gapPeakHi: Double? =
if (params.gapPeakMarginStage2 > 0.0) T + params.gapPeakMarginStage2 else null
spans = SpanUtils.mergeBimodal(
spans, sg, otsuThr, params.postMaxIdx, gapPeakHi
)
}
if (spans.isEmpty()) {
return Result(
mode, raw, sg, cfarThr, T, lowMask, rawSpans, spans,
ant = null, post = null, lowStart = null, lowEnd = null,
lowMean = null, peakMin = null, urineLen = null,
failReason = "no low-echo span"
)
}
val first = spans[0]
val s = first.start
val e = first.end
var lowSum = 0.0
for (i in s..e) lowSum += sg[i]
val lowMean = lowSum / (e - s + 1)
// py2 method_b — peak_min must satisfy BOTH (a) low_mean + margin,
// (b) >= low_echo_amp (so that wall peaks aren't picked from below T).
val peakMin = maxOf(lowMean + params.minPeakMargin, T)
var ant = WallSelect.selectWallByProminence(
sg, edge = s, searchWin = params.peakSearchWin,
peakMin = peakMin, side = WallSelect.Side.ANT, otherEdge = e,
maxCandidates = params.maxPeakCandidatesAnt,
edgeDistDecay = params.edgeDistDecay,
valleyStopRise = params.valleyStopRise,
)
// V4.1 ant — Kremkau 2017 half-amplitude rule as the PRIMARY
// ant selector. ant = lumen_start − 1 is the last sample where
// the envelope exceeds the detection threshold before
// transitioning to the hypoechoic baseline. This gives
// anatomically-consistent ant positions across clean separable
// walls, off-axis off-bladder echoes, phantom-on-floor, and
// thin-wall body phantoms where the real wall merges into the
// ringdown tail. Prominence-based ant is retained as fallback
// when (a) the lumen-edge sample is below peakMin, or (b)
// lumen_start − 1 falls inside the Fresnel near-field cutoff.
if (mode is Mode.Adaptive && params.antMinIdx > 0) {
val fb = maxOf(0, s - 1)
if (sgForMask[fb] >= peakMin && fb >= params.antMinIdx) {
ant = fb
}
}
var post = WallSelect.selectWallByProminence(
sg, edge = e, searchWin = params.peakSearchWin,
peakMin = peakMin, side = WallSelect.Side.POST, otherEdge = s,
maxCandidates = params.maxPeakCandidatesPost,
edgeDistDecay = params.edgeDistDecay,
valleyStopRise = params.valleyStopRise,
)
if (ant == null || post == null) {
return Result(
mode, raw, sg, cfarThr, T, lowMask, rawSpans, spans,
ant = null, post = null, lowStart = null, lowEnd = null,
lowMean = lowMean, peakMin = peakMin, urineLen = null,
failReason = "wall peak not found"
)
}
if (post > params.postMaxIdx) {
val backHalf = ((s + e) / 2)
val post2 = WallSelect.selectWallByProminence(
sg, edge = backHalf,
searchWin = params.postMaxIdx - backHalf,
peakMin = peakMin, side = WallSelect.Side.POST,
otherEdge = null,
maxCandidates = params.maxPeakCandidatesPost,
edgeDistDecay = params.edgeDistDecay,
valleyStopRise = params.valleyStopRise,
)
if (post2 == null) {
return Result(
mode, raw, sg, cfarThr, T, lowMask, rawSpans, spans,
ant = null, post = null, lowStart = null, lowEnd = null,
lowMean = lowMean, peakMin = peakMin, urineLen = null,
failReason = "post>POST_MAX_IDX retry failed"
)
}
post = post2
}
val urineLen = post - ant - 1
if (urineLen < params.minUrineLen) {
return Result(
mode, raw, sg, cfarThr, T, lowMask, rawSpans, spans,
ant = null, post = null, lowStart = null, lowEnd = null,
lowMean = lowMean, peakMin = peakMin, urineLen = urineLen,
failReason = "urine_len=$urineLen < 3"
)
}
return Result(
mode, raw, sg, cfarThr, T, lowMask, rawSpans, spans,
ant = ant, post = post, lowStart = s, lowEnd = e,
lowMean = lowMean, peakMin = peakMin, urineLen = urineLen,
failReason = null
)
}
}
@@ -0,0 +1,125 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Subset of study/wall_detect_verify/js/algo/geometry_v2_30deg.js (sampleToMm only).
* Full geometry (wallIdxToXyz, 12 wall points) lands in PR-6 (sphere fit).
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdConfig
import com.medithings.vesiscan.walldetect.core.WdProbe
import kotlin.math.cos
import kotlin.math.sin
object Geometry {
private const val DEG = Math.PI / 180.0
/**
* Sample index → mm depth from probe surface.
* mm = idx * dps + delayMm
*
* dps default = 1.9309 mm/sample (TB370FU 6-channel, c_eff = 1544.7 m/s, fs = 200 kHz).
* delayMm default = 6.85 mm (acoustic delay through housing + skin coupling).
*/
fun sampleToMm(
idx: Double,
dps: Double = WdConfig.DPS_DEFAULT,
delayMm: Double = WdConfig.DELAY_MM_DEFAULT
): Double = idx * dps + delayMm
/**
* V2 30° probe geometry (py4 mirror — SI-only by default).
* p = (sensor_x + d sin θ_si, 0, sensor_z + d cos θ_si)
* SI+LR (when useLr=true):
* b = (cos θ_lr · sin θ_si, sin θ_lr, cos θ_lr · cos θ_si)
* p = sensor + d · b
*/
fun wallIdxToXyz(
channel: Int,
idx: Int,
dps: Double = WdConfig.DPS_DEFAULT,
delayMm: Double = WdConfig.DELAY_MM_DEFAULT,
useLr: Boolean = WdConfig.USE_LR_TILT
): DoubleArray {
val tSi = WdProbe.DEGREE[channel] * DEG
val dist = sampleToMm(idx.toDouble(), dps, delayMm)
return if (useLr) {
val tLr = WdProbe.DEGREE_LR[channel] * DEG
val bx = cos(tLr) * sin(tSi)
val by = sin(tLr)
val bz = cos(tLr) * cos(tSi)
doubleArrayOf(
WdProbe.SENSOR_X[channel] + dist * bx,
dist * by,
WdProbe.SENSOR_Z[channel] + dist * bz
)
} else {
doubleArrayOf(
WdProbe.SENSOR_X[channel] + dist * sin(tSi),
0.0,
WdProbe.SENSOR_Z[channel] + dist * cos(tSi)
)
}
}
/**
* Visual-coordinate variant (amode_simulator convention):
* x = lateral, y = depth (forward into body, +), z = vertical
* beam dir = (sin lr · cos si, cos lr · cos si, sin si)
* SI angle is negative in PROBE.DEGREE for down-tilt → sin si < 0 → beam descends in z.
*/
fun wallIdxToXyzVisual(
channel: Int,
idx: Int,
dps: Double = WdConfig.DPS_DEFAULT,
delayMm: Double = WdConfig.DELAY_MM_DEFAULT,
useLr: Boolean = WdConfig.USE_LR_TILT
): DoubleArray {
val tSi = WdProbe.DEGREE[channel] * DEG
val tLr = if (useLr) WdProbe.DEGREE_LR[channel] * DEG else 0.0
val dist = idx * dps + delayMm
return doubleArrayOf(
WdProbe.SENSOR_X[channel] + dist * sin(tLr) * cos(tSi),
dist * cos(tLr) * cos(tSi),
WdProbe.SENSOR_Z[channel] + dist * sin(tSi)
)
}
/** A wall point: (channel, kind, xyz). kind ∈ {"ant", "post"}. */
data class WallPoint(val ch: Int, val kind: String, val xyz: DoubleArray)
/** Per-channel detection input — only ant/post indices matter. */
data class Detection(val ant: Int?, val post: Int?)
fun buildWallPoints(
detResults: List<Detection>,
dps: Double = WdConfig.DPS_DEFAULT,
delayMm: Double = WdConfig.DELAY_MM_DEFAULT,
useLr: Boolean = WdConfig.USE_LR_TILT
): List<WallPoint> {
val pts = mutableListOf<WallPoint>()
detResults.forEachIndexed { ch, r ->
r.ant?.let { pts += WallPoint(ch, "ant", wallIdxToXyz(ch, it, dps, delayMm, useLr)) }
r.post?.let { pts += WallPoint(ch, "post", wallIdxToXyz(ch, it, dps, delayMm, useLr)) }
}
return pts
}
fun buildWallPointsVisual(
detResults: List<Detection>,
dps: Double = WdConfig.DPS_DEFAULT,
delayMm: Double = WdConfig.DELAY_MM_DEFAULT,
useLr: Boolean = WdConfig.USE_LR_TILT
): List<WallPoint> {
val pts = mutableListOf<WallPoint>()
detResults.forEachIndexed { ch, r ->
r.ant?.let { pts += WallPoint(ch, "ant", wallIdxToXyzVisual(ch, it, dps, delayMm, useLr)) }
r.post?.let { pts += WallPoint(ch, "post", wallIdxToXyzVisual(ch, it, dps, delayMm, useLr)) }
}
return pts
}
}
@@ -0,0 +1,172 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/impulse_reject.js (1:1).
*
* Hampel-based lumen-aware two-pass impulse rejection — V4.1 ONLY.
*
* References
* Hampel FR. "The Influence Curve and Its Role in Robust Estimation."
* J Am Stat Assoc 69(346):383-393, 1974.
* Pearson RK, Neuvo Y, Astola J, Gabbouj M. "Generalized Hampel
* Filters." EURASIP J Adv Signal Process 2016:87, 2016.
* doi:10.1186/s13634-016-0383-6
*
* Two-pass detection workflow
* 1. Pass 1 — DetectLumenFirst (adaptive) on raw envelope → coarse (ant, post).
* 2. Range-restricted Hampel: apply ONLY to samples in [coarse_ant+1,
* coarse_post-1]. Outliers replaced by local median. Walls untouched.
* 3. Pass 2 — DetectLumenFirst on cleaned envelope → final (ant, post).
*
* The two-pass structure is the architectural strength of the V4.1
* pipeline: lumen-bump cleaning needs to know where the lumen is; without
* Pass 1 wall localisation the Hampel filter has no anchor.
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
object ImpulseReject {
const val DEFAULT_WIN = 11
const val DEFAULT_K = 3.0
data class FlaggedSample(
val i: Int,
val original: Double,
val replaced: Double,
val deviation: Double,
val sigma: Double
)
data class HampelOutput(val out: DoubleArray, val flagged: List<FlaggedSample>)
private fun localMedian(x: DoubleArray, lo: Int, hi: Int): Pair<Double, Double> {
val w = DoubleArray(hi - lo + 1) { x[lo + it] }
val s = w.sortedArray()
val med = s[s.size / 2]
val dev = DoubleArray(s.size) { abs(w[it] - med) }
val ds = dev.sortedArray()
val mad = ds[ds.size / 2]
return Pair(med, mad)
}
/** Plain Hampel filter — full envelope. */
fun hampel(x: DoubleArray, win: Int = DEFAULT_WIN, k: Double = DEFAULT_K): HampelOutput {
val n = x.size
val half = (win - 1) / 2
val out = DoubleArray(n)
val flagged = mutableListOf<FlaggedSample>()
for (i in 0 until n) {
val lo = maxOf(0, i - half)
val hi = minOf(n - 1, i + half)
val (med, mad) = localMedian(x, lo, hi)
val sigma = 1.4826 * mad
if (sigma > 0.0 && abs(x[i] - med) > k * sigma) {
out[i] = med
flagged.add(FlaggedSample(i, x[i], med, x[i] - med, sigma))
} else {
out[i] = x[i]
}
}
return HampelOutput(out, flagged)
}
/**
* Range-restricted Hampel: only operates on samples i ∈ [lo, hi].
* Samples outside copy through unchanged — wall preservation guarantee.
*/
fun hampelInRange(
x: DoubleArray,
lo: Int,
hi: Int,
win: Int = DEFAULT_WIN,
k: Double = DEFAULT_K
): HampelOutput {
val out = x.copyOf()
val flagged = mutableListOf<FlaggedSample>()
if (lo >= hi) return HampelOutput(out, flagged)
val half = (win - 1) / 2
for (i in lo..hi) {
val wlo = maxOf(0, i - half)
val whi = minOf(x.size - 1, i + half)
val (med, mad) = localMedian(x, wlo, whi)
val sigma = 1.4826 * mad
if (sigma > 0.0 && abs(x[i] - med) > k * sigma) {
out[i] = med
flagged.add(FlaggedSample(i, x[i], med, x[i] - med, sigma))
}
}
return HampelOutput(out, flagged)
}
data class TwoPassResult(
val coarse: DetectLumenFirst.Result,
val refined: DetectLumenFirst.Result,
val cleanedEnvelope: DoubleArray,
val bumpsRemoved: List<Int>,
val flagged: List<FlaggedSample>
)
/**
* Lumen-aware two-pass V4.1 detection.
* - Pass 1: standard adaptive detect on raw envelope.
* - Hampel within [coarse_ant + 1 + ⌊W/2⌋, coarse_post − 1 − ⌊W/2⌋].
* - Pass 2: re-detect on cleaned envelope.
*
* §3.6 Symmetric edge-bias buffer (Pearson–Neuvo 2016 §4.2):
* the Hampel local window of width W has a step-discontinuity ripple
* region exactly ⌊W/2⌋ samples wide on each side of a wall transition;
* shrinking the operating range by ⌊W/2⌋ from each lumen boundary
* guarantees the local window never straddles a wall sample, so the
* MAD does not inflate and lumen-edge samples are not falsely flagged.
* The buffer width is the closed-form derivative of the existing W
* parameter — no new constants. See PUBLICATION-ROADMAP.md §A.4 for
* the manuscript-side framing.
*/
fun detectWithLumenClean(
rawAdc: DoubleArray,
mode: DetectLumenFirst.Mode = DetectLumenFirst.Mode.Adaptive,
win: Int = DEFAULT_WIN,
k: Double = DEFAULT_K,
params: DetectLumenFirst.ParamSet = DetectLumenFirst.V41_PARAMS
): TwoPassResult {
val coarse = DetectLumenFirst.detect(rawAdc, mode, params)
if (coarse.ant == null || coarse.post == null) {
return TwoPassResult(coarse, coarse, rawAdc.copyOf(), emptyList(), emptyList())
}
// §3.6 — symmetric edge-bias buffer, closed-form from win.
val half = (win - 1) / 2
val lo = coarse.ant + 1 + half
val hi = coarse.post - 1 - half
if (lo >= hi) {
// Lumen too narrow for any safe Hampel window — pass through.
return TwoPassResult(coarse, coarse, rawAdc.copyOf(), emptyList(), emptyList())
}
val cleaned = hampelInRange(rawAdc, lo, hi, win, k)
if (cleaned.flagged.isEmpty()) {
return TwoPassResult(coarse, coarse, rawAdc.copyOf(), emptyList(), emptyList())
}
val refined = DetectLumenFirst.detect(cleaned.out, mode, params)
return TwoPassResult(
coarse = coarse,
refined = refined,
cleanedEnvelope = cleaned.out,
bumpsRemoved = cleaned.flagged.map { it.i },
flagged = cleaned.flagged
)
}
fun detectWithLumenClean(
rawAdc: IntArray,
mode: DetectLumenFirst.Mode = DetectLumenFirst.Mode.Adaptive,
win: Int = DEFAULT_WIN,
k: Double = DEFAULT_K,
params: DetectLumenFirst.ParamSet = DetectLumenFirst.V41_PARAMS
): TwoPassResult {
val raw = DoubleArray(rawAdc.size) { rawAdc[it].toDouble() }
return detectWithLumenClean(raw, mode, win, k, params)
}
}
@@ -0,0 +1,121 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/median_filter.js (1:1).
*
* 1-D running-median pre-filter for the V4.1 adaptive path.
*
* Why this exists
* On the 150 mL Japan body phantom (CCK 2026-04-28) channels CH0/CH1
* exhibited intermittent detection failure. Diagnosis showed isolated
* high-amplitude speckle bumps INSIDE the lumen at alternating samples
* (37, 39, 41, 43, 45 in CH0). Five bumps interspersed with a
* hypoechoic baseline (~830 ADC) just above the OS-CFAR threshold
* (T ≈ 887 ADC):
*
* window-of-11 around bump 1216:
* sorted = {850, 850, 850, 850, 850, 850, 979, 992, 1082, 1142, 1216}
* median = 850 ← the 6th value
* MAD = median{|x − 850|} = 0 ← the 6th deviation is zero
* σ = 1.4826 · MAD = 0
*
* Hampel's outlier test reads
* if (sigma > 0 && |x[i] − med| > k · sigma) flag x[i];
* With σ = 0 the guard short-circuits and **no bump is flagged**.
* This is the degenerate-MAD masking failure described in
*
* Davies, L. and Gather, U. "The identification of multiple
* outliers." J Am Stat Assoc 88(423):782–792, 1993.
*
* When the outlier density inside the filter window exceeds ~50 %,
* the median is pulled into the bump cluster and MAD collapses,
* making clustered impulses invisible to Hampel.
*
* Method (literature)
* Tukey, J.W. "Nonlinear (nonsuperposable) methods for smoothing data."
* Cong Rec 1974 EASCON, p673. Original running-median proposal.
* Justusson, B.I. "Median filtering: Statistical properties." In:
* Two-Dimensional Digital Signal Processing II (Topics in Applied
* Physics 43), Springer 1981, p161–196. Convergence and root-signal
* theory: impulses ≤ ⌊win/2⌋ samples wide are guaranteed removed.
* Loizou, C.P. and Pattichis, C.S. "Despeckle Filtering Algorithms and
* Software for Ultrasound Imaging." Synthesis Lectures on Algorithms
* and Software in Engineering, Morgan & Claypool 2008. Median is the
* reference baseline despeckle method against which adaptive filters
* (Lee 1980, Frost 1982) are compared.
*
* Width selection (7)
* For a length-W running median, isolated impulses up to ⌊W/2⌋
* consecutive samples are absorbed (Justusson 1981 §2). W = 7 absorbs
* 1–3-sample bumps — exactly matching the Burckhardt 1978 prediction
* of 1–3-sample speckle peaks in hypoechoic regions, and complementary
* to MERGE_GAP_MAX = 5 / GAP_PEAK_MARGIN = 50 (which handles wider
* low-amplitude speckle clusters). Real bladder-wall echoes span
* ≥ 4 samples and pass through the filter unchanged: by Justusson 1981
* Theorem 2.3 every plateau of length ≥ ⌈W/2⌉ + 1 = 4 is a fixed point
* of the W = 7 median.
*
* Edge handling
* Symmetric reflection (Gonzalez-Woods 2017 §3.4) at both ends.
*/
package com.medithings.vesiscan.walldetect.algo
object MedianFilter {
fun runningMedian(x: DoubleArray, win: Int): DoubleArray {
if (win < 2) return x.copyOf()
val half = (win - 1) / 2
val n = x.size
val out = DoubleArray(n)
val buf = DoubleArray(win)
fun sample(i: Int): Double {
var k = i
if (k < 0) k = -k - 1
if (k >= n) k = 2 * n - k - 1
if (k < 0) k = 0
if (k >= n) k = n - 1
return x[k]
}
for (i in 0 until n) {
for (j in 0 until win) buf[j] = sample(i + j - half)
// Insertion sort — win is typically 7.
for (a in 1 until win) {
val v = buf[a]
var b = a - 1
while (b >= 0 && buf[b] > v) {
buf[b + 1] = buf[b]
b--
}
buf[b + 1] = v
}
out[i] = buf[half]
}
return out
}
/**
* Iterated running-median → fixed-point (root signal).
* Justusson 1981 §4: a finite number of passes drives any input to
* a fixed point of the W-median operator. For W = 7 on the 150 mL
* Japan body phantom CH0 trace (5 alternating bumps spanning 9
* samples) convergence is reached in two passes; we cap at 4 for
* safety.
*/
fun runningMedianRoot(x: DoubleArray, win: Int, maxIters: Int = 4): DoubleArray {
var cur = runningMedian(x, win)
for (it in 1 until maxIters) {
val next = runningMedian(cur, win)
var same = true
for (i in cur.indices) {
if (cur[i] != next[i]) { same = false; break }
}
cur = next
if (same) break
}
return cur
}
}
@@ -0,0 +1,66 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/morph_close.js (1:1).
*
* 1-D morphological closing on the binary lumen mask.
*
* References
* Serra J. "Image Analysis and Mathematical Morphology." Academic
* Press, 1982.
* Soille P. "Morphological Image Analysis: Principles and Applications."
* 2nd ed., Springer, 2003.
*
* Closing(M, B) = Erode(Dilate(M, B), B)
*
* NOTE: Empirical testing on the bench's 3-capture set showed that a
* length-5 structuring element causes regressions on corner / center
* captures (CH3 fails to detect because closing merges across real wall
* transitions). The amplitude-aware MERGE_GAP_MAX=5 + GAP_PEAK_MARGIN=50
* tweak in DetectLumenFirst.kt is the chosen V4.1 default. This module is
* kept available for future use (e.g. larger structuring elements with
* length-aware decimation).
*/
package com.medithings.vesiscan.walldetect.algo
object MorphClose {
private fun halfBefore(n: Int) = (n - 1) / 2
private fun halfAfter(n: Int) = n / 2
/** Dilation: TRUE if any sample within the structuring window is TRUE. */
fun dilate(mask: BooleanArray, n: Int): BooleanArray {
val N = mask.size
val hb = halfBefore(n); val ha = halfAfter(n)
val out = BooleanArray(N)
for (i in 0 until N) {
val lo = maxOf(0, i - hb); val hi = minOf(N - 1, i + ha)
for (k in lo..hi) {
if (mask[k]) { out[i] = true; break }
}
}
return out
}
/** Erosion: TRUE only if every sample within the structuring window is TRUE. */
fun erode(mask: BooleanArray, n: Int): BooleanArray {
val N = mask.size
val hb = halfBefore(n); val ha = halfAfter(n)
val out = BooleanArray(N) { true }
for (i in 0 until N) {
val lo = maxOf(0, i - hb); val hi = minOf(N - 1, i + ha)
for (k in lo..hi) {
if (!mask[k]) { out[i] = false; break }
}
}
return out
}
/** Closing = dilation ∘ erosion. Fills gaps shorter than n. */
fun closeMask(mask: BooleanArray, n: Int): BooleanArray = erode(dilate(mask, n), n)
/** Opening = erosion ∘ dilation. Removes islands shorter than n. */
fun openMask(mask: BooleanArray, n: Int): BooleanArray = dilate(erode(mask, n), n)
}
@@ -0,0 +1,111 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/multi_mode_bv.js (1:1).
*
* Confidence-gated mode selection (§11):
* Mode A : multi-channel Kasa→LM fit, when N_trust ≥ N_min
* Mode B : single-channel chord-based BV, median over trusted channels
* (1 ≤ N_trust < N_min)
* Mode C : no measurement (N_trust = 0)
*
* Trust criterion: s_contrast(channel) ≥ τ (default 80, "low tier and above").
* N_min = 3 channels.
*
* NOTE: trust here is driven by ContrastAux (legacy s_contrast), NOT BModeScore.
* Two trust criteria coexist — UI tier comes from BModeScore, mode selection
* comes from ContrastAux. See DESIGN §1 footnote.
*/
package com.medithings.vesiscan.walldetect.algo
object MultiModeBv {
/** Per-channel detection input — needs sg + ant + post for s_contrast. */
data class Detection(val sg: DoubleArray, val ant: Int?, val post: Int?)
/** Per-channel s_contrast / tier / trust outcome. */
data class PerChannel(
val ch: Int,
val ant: Int?,
val post: Int?,
val contrast: Double?,
val tier: String, // "high" | "moderate" | "low" | "zero" | "fail" | "—"
val trust: Boolean
)
/** selectMode return type. */
data class Result(
val mode: String, // "A" | "B" | "C"
val bv: Double?,
val r: Double?,
val nTrust: Int,
val perChannel: List<PerChannel>,
val trustedSet: Map<Int, Boolean>,
val fit: SphereFit2Step.FitResult? = null,
val agg: SingleChannelBv.AggResult? = null,
val reason: String? = null
)
fun selectMode(
detections: List<Detection>,
tau: Double = 80.0,
nMin: Int = 3,
useLr: Boolean = false,
reduction: SingleChannelBv.Reduction = SingleChannelBv.Reduction.MEDIAN
): Result {
// Step 1 — compute s_contrast per channel
val perCh: List<PerChannel> = detections.mapIndexed { ch, r ->
if (r.ant == null || r.post == null) {
PerChannel(ch, null, null, null, "fail", false)
} else {
val a = ContrastAux.farPostContrast(r.sg, r.ant, r.post)
val tier = ContrastAux.classify(a?.contrast).tier
val trust = a != null && a.contrast >= tau
PerChannel(ch, r.ant, r.post, a?.contrast, tier, trust)
}
}
val trustedSet: Map<Int, Boolean> = perCh.associate { it.ch to it.trust }
val nTrust = perCh.count { it.trust }
// Step 2 — mode selection
if (nTrust == 0) {
return Result(
mode = "C", bv = null, r = null, nTrust = 0,
perChannel = perCh, trustedSet = trustedSet,
reason = "no channel passes confidence floor"
)
}
if (nTrust >= nMin) {
// Mode A: multi-channel Kasa→LM on visual coords
val trustedDet = detections.mapIndexed { ch, det ->
if (trustedSet[ch] == true) Geometry.Detection(det.ant, det.post)
else Geometry.Detection(null, null)
}
val wp = Geometry.buildWallPointsVisual(trustedDet, useLr = useLr)
val fit = SphereFit2Step.fit2Step(wp.map { it.xyz }, SphereFit2Step.Mode.AUTO)
?: return Result(
mode = "C", bv = null, r = null, nTrust = nTrust,
perChannel = perCh, trustedSet = trustedSet,
reason = "sphere fit failed despite N_trust ≥ N_min"
)
val rMm = fit.lmR
val bv = BvFromSphere.bvFromSphere(rMm)
return Result(
mode = "A", bv = bv, r = rMm, nTrust = nTrust,
perChannel = perCh, trustedSet = trustedSet, fit = fit
)
}
// Mode B: single-channel chord aggregation
val perBV = SingleChannelBv.chordBVPerChannel(
detections.map { Geometry.Detection(it.ant, it.post) }
)
val agg = SingleChannelBv.aggregate(perBV, trustedSet, reduction)
return Result(
mode = "B", bv = agg.bv, r = agg.r, nTrust = nTrust,
perChannel = perCh, trustedSet = trustedSet, agg = agg
)
}
}
@@ -0,0 +1,181 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Direct port of py2/for_app_share/span_utils.py otsu_1d().
* 1-D Otsu binary threshold — k-means(k=2) equivalent on a histogram.
* Used by detect_low_echo to derive an adaptive LOW_ECHO_AMP per signal.
*/
package com.medithings.vesiscan.walldetect.algo
object Otsu {
/** Otsu threshold + separability — Python `otsu_1d(values)` 의 (threshold, sep) tuple 호환. */
data class OtsuResult(val threshold: Double, val separability: Double)
fun otsu1dWithSeparability(values: DoubleArray, nBins: Int = 64): OtsuResult {
if (values.isEmpty()) return OtsuResult(0.0, 0.0)
var lo = Double.POSITIVE_INFINITY
var hi = Double.NEGATIVE_INFINITY
var sum = 0.0
for (v in values) {
if (v < lo) lo = v
if (v > hi) hi = v
sum += v
}
if (values.size == 1 || lo == hi) return OtsuResult(sum / values.size, 0.0)
val hist = IntArray(nBins)
val width = (hi - lo) / nBins
for (v in values) {
var b = ((v - lo) / width).toInt()
if (b >= nBins) b = nBins - 1
if (b < 0) b = 0
hist[b]++
}
val total = values.size
val centers = DoubleArray(nBins) { lo + (it + 0.5) * width }
val p = DoubleArray(nBins) { hist[it].toDouble() / total }
var muT = 0.0
for (i in 0 until nBins) muT += p[i] * centers[i]
var sigmaT = 0.0
for (i in 0 until nBins) {
val d = centers[i] - muT
sigmaT += p[i] * d * d
}
if (sigmaT <= 1e-12) return OtsuResult(muT, 0.0)
var cumP = 0.0
var cumMP = 0.0
var bestT = 0
var bestSigmaB = Double.NEGATIVE_INFINITY
for (t in 0 until nBins - 1) {
cumP += p[t]
cumMP += p[t] * centers[t]
val w0 = cumP
val w1 = 1.0 - w0
if (w0 <= 1e-6 || w1 <= 1e-6) continue
val m0 = cumMP / w0
val m1 = (muT - cumMP) / w1
val sigmaB = w0 * w1 * (m0 - m1) * (m0 - m1)
if (sigmaB > bestSigmaB) {
bestSigmaB = sigmaB
bestT = t
}
}
val sep = (bestSigmaB / sigmaT).coerceIn(0.0, 1.0)
return OtsuResult(centers[bestT], sep)
}
/**
* 1-D Otsu threshold over `values`. Returns 0.0 for empty input,
* `values.mean()` for single-value or constant input. n_bins default 64
* (32-128 range is stable per py2 docstring).
*/
fun otsu1d(values: DoubleArray, nBins: Int = 64): Double {
if (values.isEmpty()) return 0.0
var lo = Double.POSITIVE_INFINITY
var hi = Double.NEGATIVE_INFINITY
var sum = 0.0
for (v in values) {
if (v < lo) lo = v
if (v > hi) hi = v
sum += v
}
if (values.size == 1 || lo == hi) return sum / values.size
// Histogram
val hist = IntArray(nBins)
val width = (hi - lo) / nBins
for (v in values) {
// numpy.histogram: rightmost bin includes the right edge
var b = ((v - lo) / width).toInt()
if (b >= nBins) b = nBins - 1
if (b < 0) b = 0
hist[b]++
}
val total = values.size
if (total == 0) return sum / values.size
val centers = DoubleArray(nBins) { lo + (it + 0.5) * width }
val p = DoubleArray(nBins) { hist[it].toDouble() / total }
var cumP = 0.0
var cumMP = 0.0
val cumPArr = DoubleArray(nBins)
val cumMPArr = DoubleArray(nBins)
for (i in 0 until nBins) {
cumP += p[i]
cumMP += p[i] * centers[i]
cumPArr[i] = cumP
cumMPArr[i] = cumMP
}
val totalM = cumMPArr[nBins - 1]
var bestT = 0
var bestSigma = Double.NEGATIVE_INFINITY
for (t in 0 until nBins - 1) {
val w0 = cumPArr[t]
val w1 = 1.0 - w0
if (w0 <= 1e-6 || w1 <= 1e-6) continue
val m0 = cumMPArr[t] / w0
val m1 = (totalM - cumMPArr[t]) / w1
val sigmaB = w0 * w1 * (m0 - m1) * (m0 - m1)
if (sigmaB > bestSigma) {
bestSigma = sigmaB
bestT = t
}
}
return centers[bestT]
}
/**
* 1-ADC-resolution Otsu (Otsu 1979) — one histogram bin per integer
* ADC value. The 64-bin variant above is fast for visualisation but
* its bin width on a ~1200-ADC envelope is ~19 ADC, which can shift
* the threshold by ±10–20 ADC vs an unbinned histogram. For
* gap-merge bimodal discrimination (V4.1 §6.5e Stage 2) that ±10
* ADC is enough to flip the decision when a wall peak sits within
* ~50 ADC of the cluster ceiling.
*/
fun otsu1dInteger(values: DoubleArray): Double {
if (values.isEmpty()) return 0.0
var lo = Double.POSITIVE_INFINITY
var hi = Double.NEGATIVE_INFINITY
for (v in values) {
if (v < lo) lo = v
if (v > hi) hi = v
}
val loInt = kotlin.math.floor(lo).toInt()
val hiInt = kotlin.math.ceil(hi).toInt()
if (hiInt == loInt) return loInt.toDouble()
val bins = hiInt - loInt + 1
val hist = IntArray(bins)
for (v in values) {
var k = (kotlin.math.round(v) - loInt).toInt()
if (k < 0) k = 0 else if (k >= bins) k = bins - 1
hist[k]++
}
val total = values.size
var sumAll = 0.0
for (i in 0 until bins) sumAll += i * hist[i]
var sumB = 0.0
var wB = 0
var maxVar = -1.0
var bestI = 0
for (i in 0 until bins) {
wB += hist[i]
if (wB == 0) continue
val wF = total - wB
if (wF == 0) break
sumB += i * hist[i]
val mB = sumB / wB
val mF = (sumAll - sumB) / wF
val v = wB.toDouble() * wF * (mB - mF) * (mB - mF)
if (v > maxVar) {
maxVar = v
bestI = i
}
}
return (loInt + bestI).toDouble()
}
}
@@ -0,0 +1,82 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/peak_detection.js (1:1).
* Strict local maxima + prominence + valley walks (py4/peak_detection style).
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.max
import kotlin.math.min
object PeakDetection {
/**
* Strict local maxima on x[lo..hi). Plateau-aware: a flat plateau is
* reported as a single peak at its midpoint (rounded down to int).
*/
fun findPeaks1D(
x: DoubleArray,
lo: Int = 0,
hi: Int = x.size,
minHeight: Double? = null
): IntArray {
val peaks = mutableListOf<Int>()
var i = lo + 1
while (i < hi - 1) {
if (x[i - 1] < x[i]) {
var j = i
while (j < hi - 1 && x[j + 1] == x[j]) j++
if (j < hi - 1 && x[j + 1] < x[j]) {
peaks += ((i + j) / 2) // integer division
}
i = j + 1
} else {
i++
}
}
return if (minHeight == null) peaks.toIntArray()
else peaks.filter { x[it] >= minHeight }.toIntArray()
}
/**
* Prominence with valley-walk and break-rise stop (py4 style on log envelope).
*/
fun prominence(x: DoubleArray, p: Int, maxDist: Int = 25, breakRise: Double = 1.0): Double {
val n = x.size
var vr = x[p]
for (k in (p + 1) until min(n, p + maxDist + 1)) {
if (x[k] < vr) vr = x[k]
else if (x[k] > vr + breakRise) break
}
var vl = x[p]
for (k in (p - 1) downTo max(-1, p - maxDist - 1) + 1) {
if (x[k] < vl) vl = x[k]
else if (x[k] > vl + breakRise) break
}
return x[p] - max(vl, vr)
}
/** Right valley walk on raw amplitude (py2 style). */
fun rightValley(sg: DoubleArray, p: Int, maxDist: Int = 20, breakRise: Double = 10.0): Double {
val n = sg.size
var v = sg[p]
for (i in (p + 1) until min(n, p + maxDist)) {
if (sg[i] < v) v = sg[i]
else if (sg[i] > v + breakRise) break
}
return v
}
/** Left valley walk on raw amplitude (py2 style). */
fun leftValley(sg: DoubleArray, p: Int, maxDist: Int = 20, breakRise: Double = 10.0): Double {
var v = sg[p]
for (i in (p - 1) downTo max(0, p - maxDist)) {
if (sg[i] < v) v = sg[i]
else if (sg[i] > v + breakRise) break
}
return v
}
}
@@ -0,0 +1,139 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Direct port of py2 find_plateau_with_peaks_v2 (study/py2/labeling.py L332).
* This replaces the JS lumen-first detector for V2 (per user, 2026-04-28).
*
* Pipeline:
* 1. Sliding plateau score (sliding_scores_1d, win=5)
* 2. plateau_mask = score ≤ quantile(score, 0.7)
* 3. spans → length filter [L_min=5, L_max=30] → merge_overlapping
* 4. Detrended signal (uniform_filter1d, size=15)
* 5. find_peaks(detrended, prominence=10, distance=3)
* 6. For each plateau, find left/right peak with peak_contrast ≥ 0.02
* 7. final segment = (left_peak+1, right_peak-1)
* 8. merge_overlapping + merge_adjacent_regions
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
import kotlin.math.max
import kotlin.math.min
object PlateauDetector {
// py2 config defaults (PLATEAU_*)
const val WIN = 5
const val SCORE_Q = 0.7
const val L_MIN = 5
const val L_MAX = 30
const val PROMINENCE = 10.0
const val DETREND_WIN = 15
const val PEAK_CONTRAST = 0.02
/** A single plateau with its left/right peak (anterior/posterior wall). */
data class WallPair(
val ant: Int, // = left_peak (sample idx)
val post: Int, // = right_peak
val plateauStart: Int, // detected plateau (inclusive)
val plateauEnd: Int,
)
data class Result(
val raw: DoubleArray,
val score: DoubleArray,
val plateauSegs: List<IntRange>, // raw plateaus (after L_min/L_max filter)
val finalSegments: List<IntRange>, // post peak validation + merging
val wallPairs: List<WallPair>, // primary output for V2Detector
)
fun detect(
rawAdc: DoubleArray,
win: Int = WIN,
scoreQ: Double = SCORE_Q,
lMin: Int = L_MIN,
lMax: Int = L_MAX,
prominence: Double = PROMINENCE,
detrendWin: Int = DETREND_WIN,
peakContrast: Double = PEAK_CONTRAST,
): Result {
val x = rawAdc.copyOf()
val n = x.size
// 1. Plateau detection
val sliding = Py2Helpers.slidingScores1d(x, win)
val score = sliding.score
val thr = Py2Helpers.npQuantile(score, scoreQ)
val plateauMask = BooleanArray(n) { score[it] <= thr }
val raw = Py2Helpers.maskToSegments(plateauMask)
val filtered = raw.filter { (it.last - it.first + 1) in lMin..lMax }
val plateauSegs = Py2Helpers.mergeOverlappingSegments(filtered)
if (plateauSegs.isEmpty()) {
return Result(x, score, emptyList(), emptyList(), emptyList())
}
// 2. Detrend
val trend = Py2Helpers.uniformFilter1d(x, detrendWin)
val xDetrended = DoubleArray(n) { x[it] - trend[it] }
val allPeaks = Py2Helpers.findPeaks(xDetrended, prominence, distance = 3)
// 3. Per-plateau peak validation
val finalSegs = mutableListOf<IntRange>()
val pairs = mutableListOf<WallPair>()
for (plateau in plateauSegs) {
val ps = plateau.first
val pe = plateau.last
var sum = 0.0
for (k in ps..pe) sum += x[k]
val plateauMean = sum / (pe - ps + 1)
val searchRange = max(15, pe - ps)
// Left peak — closest to plateau, scanning right→left
val leftLo = max(0, ps - searchRange)
val leftCands = allPeaks.filter { it in leftLo until ps }
var leftPeak: Int? = null
for (lp in leftCands.reversed()) {
if ((x[lp] - plateauMean) / (abs(plateauMean) + 1e-12) >= peakContrast) {
leftPeak = lp
break
}
}
// Right peak — closest to plateau, scanning left→right
val rightHi = min(n, pe + searchRange + 1)
val rightCands = allPeaks.filter { it in (pe + 1) until rightHi }
var rightPeak: Int? = null
for (rp in rightCands) {
if ((x[rp] - plateauMean) / (abs(plateauMean) + 1e-12) >= peakContrast) {
rightPeak = rp
break
}
}
if (leftPeak == null || rightPeak == null) continue
val newS = leftPeak + 1
val newE = rightPeak - 1
if (newE <= newS) continue
finalSegs += newS..newE
pairs += WallPair(ant = leftPeak, post = rightPeak, plateauStart = ps, plateauEnd = pe)
}
var merged = Py2Helpers.mergeOverlappingSegments(finalSegs)
merged = Py2Helpers.mergeAdjacentRegions(x, merged)
return Result(
raw = x,
score = score,
plateauSegs = plateauSegs,
finalSegments = merged,
wallPairs = pairs,
)
}
}
@@ -0,0 +1,302 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Direct port of py2 helpers used by find_plateau_with_peaks_v2:
* - robust_z (labeling.py L26)
* - sliding_scores_1d (labeling.py L34)
* - mask_to_segments (labeling.py L68)
* - merge_overlapping_segments (labeling.py L87)
* - merge_adjacent_regions (labeling.py L101)
* - uniform_filter1d (scipy.ndimage equivalent — mode='nearest')
* - find_peaks (scipy.signal equivalent — prominence + distance)
*
* NumPy median (even-length → mean of two middle values) is reproduced
* faithfully, since robust_z drives the per-sample plateau score.
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
import kotlin.math.ceil
import kotlin.math.floor
import kotlin.math.max
import kotlin.math.min
import kotlin.math.sqrt
object Py2Helpers {
// ─────────────────────────────────────────────────────────
// numpy.median (even N → mean of two middle values)
// ─────────────────────────────────────────────────────────
fun npMedian(values: DoubleArray): Double {
if (values.isEmpty()) return Double.NaN
val s = values.copyOf().also { it.sort() }
val n = s.size
return if (n % 2 == 0) (s[n / 2 - 1] + s[n / 2]) / 2.0 else s[n / 2]
}
/** Numpy linear-interpolation quantile. */
fun npQuantile(values: DoubleArray, q: Double): Double {
if (values.isEmpty()) return Double.NaN
val s = values.copyOf().also { it.sort() }
val pos = q * (s.size - 1)
val lo = floor(pos).toInt()
val hi = ceil(pos).toInt()
if (lo == hi) return s[lo]
return s[lo] + (s[hi] - s[lo]) * (pos - lo)
}
// ─────────────────────────────────────────────────────────
// robust_z — (a − median) / (1.4826 · MAD)
// ─────────────────────────────────────────────────────────
fun robustZ(a: DoubleArray): DoubleArray {
val med = npMedian(a)
val abs_dev = DoubleArray(a.size) { abs(a[it] - med) }
val mad = npMedian(abs_dev) + 1e-12
return DoubleArray(a.size) { (a[it] - med) / (1.4826 * mad) }
}
// ─────────────────────────────────────────────────────────
// sliding_scores_1d
// ─────────────────────────────────────────────────────────
data class SlidingScores(
val score: DoubleArray,
val flat: DoubleArray,
val slope: DoubleArray,
val low: DoubleArray,
)
fun slidingScores1d(x: DoubleArray, winInput: Int): SlidingScores {
require(winInput >= 3) { "win >= 3 required" }
var win = winInput
if (win % 2 == 0) win += 1
val half = win / 2
val t = x.size
// Edge padding (numpy: mode='edge')
val xp = DoubleArray(t + 2 * half)
for (i in 0 until half) xp[i] = x[0]
for (i in 0 until t) xp[half + i] = x[i]
for (i in 0 until half) xp[half + t + i] = x[t - 1]
// tt = arange(win) − mean
val tt = DoubleArray(win) { it.toDouble() }
val meanTt = tt.average()
for (i in tt.indices) tt[i] -= meanTt
var denom = 0.0
for (v in tt) denom += v * v
denom += 1e-12
val flat = DoubleArray(t)
val slope = DoubleArray(t)
val low = DoubleArray(t)
val w = DoubleArray(win)
for (i in 0 until t) {
for (k in 0 until win) w[k] = xp[i + k]
val mu = w.average()
var sq = 0.0
for (v in w) { val d = v - mu; sq += d * d }
flat[i] = sqrt(sq / win)
low[i] = npQuantile(w, 0.2)
var num = 0.0
for (k in 0 until win) num += tt[k] * (w[k] - mu)
slope[i] = abs(num / denom)
}
val rzLow = robustZ(low)
val rzFlat = robustZ(flat)
val rzSlope = robustZ(slope)
val score = DoubleArray(t) { rzLow[it] + rzFlat[it] + rzSlope[it] }
return SlidingScores(score, flat, slope, low)
}
// ─────────────────────────────────────────────────────────
// mask_to_segments / merge_overlapping_segments
// ─────────────────────────────────────────────────────────
fun maskToSegments(mask: BooleanArray): List<IntRange> {
val segs = mutableListOf<IntRange>()
var inRun = false
var s = 0
for (i in mask.indices) {
if (mask[i] && !inRun) { inRun = true; s = i }
else if (!mask[i] && inRun) { inRun = false; segs += s..(i - 1) }
}
if (inRun) segs += s..(mask.size - 1)
return segs
}
fun mergeOverlappingSegments(segments: List<IntRange>): List<IntRange> {
if (segments.isEmpty()) return emptyList()
val sorted = segments.sortedBy { it.first }
val merged = mutableListOf(sorted[0])
for (seg in sorted.drop(1)) {
val last = merged.last()
if (seg.first <= last.last + 1) {
merged[merged.lastIndex] = last.first..max(last.last, seg.last)
} else {
merged.add(seg)
}
}
return merged
}
// ─────────────────────────────────────────────────────────
// uniform_filter1d (scipy.ndimage; mode='nearest' = edge padding)
// ─────────────────────────────────────────────────────────
fun uniformFilter1d(x: DoubleArray, size: Int): DoubleArray {
require(size >= 1)
val t = x.size
val half = size / 2
// scipy uniform_filter1d centers the window; for odd size half=size/2.
// For even size, the bias is handled differently — we follow scipy's
// "origin=0" default which effectively uses window [i-half, i+half-(1 if even else 0)].
val left = half
val right = size - 1 - left
val xp = DoubleArray(t + left + right)
for (i in 0 until left) xp[i] = x[0]
for (i in 0 until t) xp[left + i] = x[i]
for (i in 0 until right) xp[left + t + i] = x[t - 1]
val out = DoubleArray(t)
for (i in 0 until t) {
var s = 0.0
for (k in 0 until size) s += xp[i + k]
out[i] = s / size
}
return out
}
// ─────────────────────────────────────────────────────────
// find_peaks (scipy.signal — prominence + distance)
//
// Implementation:
// 1. Strict local maxima with plateau handling (midpoint).
// 2. Per-peak prominence: walk left/right until hitting a higher sample
// or array edge, take the minimum encountered; prominence = peak −
// max(left_min, right_min).
// 3. Drop peaks below `prominence`.
// 4. Apply `distance`: greedy keep, prefer larger prominence.
// ─────────────────────────────────────────────────────────
fun findPeaks(x: DoubleArray, prominence: Double, distance: Int): IntArray {
val t = x.size
if (t < 3) return IntArray(0)
// 1. Local maxima
val candidates = mutableListOf<Int>()
var i = 1
while (i < t - 1) {
if (x[i - 1] < x[i]) {
var j = i
while (j < t - 1 && x[j + 1] == x[j]) j++
if (j < t - 1 && x[j + 1] < x[j]) {
candidates.add((i + j) / 2)
}
i = j + 1
} else {
i++
}
}
if (candidates.isEmpty()) return IntArray(0)
// 2. Prominence
val proms = DoubleArray(candidates.size)
for ((idx, p) in candidates.withIndex()) {
var leftMin = x[p]
var k = p - 1
while (k >= 0) {
if (x[k] > x[p]) break
if (x[k] < leftMin) leftMin = x[k]
k--
}
var rightMin = x[p]
k = p + 1
while (k < t) {
if (x[k] > x[p]) break
if (x[k] < rightMin) rightMin = x[k]
k++
}
proms[idx] = x[p] - max(leftMin, rightMin)
}
// 3. Filter by prominence
val filtered = candidates.indices
.filter { proms[it] >= prominence }
.map { candidates[it] to proms[it] }
if (distance <= 0) return filtered.map { it.first }.sorted().toIntArray()
// 4. Distance constraint (greedy by prominence desc)
val byProm = filtered.sortedByDescending { it.second }
val keep = mutableListOf<Int>()
for ((idx, _) in byProm) {
if (keep.none { abs(it - idx) < distance }) keep.add(idx)
}
return keep.sorted().toIntArray()
}
// ─────────────────────────────────────────────────────────
// merge_adjacent_regions — py2 labeling.py L101
// ─────────────────────────────────────────────────────────
fun mergeAdjacentRegions(
x: DoubleArray,
segments: List<IntRange>,
maxGapLen: Int = 12,
relTol: Double = 0.05,
refWin: Int = 5,
maxMergedLen: Int = 50,
): List<IntRange> {
if (segments.size < 2) return segments
val merged = mutableListOf(segments[0])
for (seg in segments.drop(1)) {
val prev = merged.last()
val prevS = prev.first
val prevE = prev.last
val currS = seg.first
val currE = seg.last
val gapStart = prevE + 1
val gapEnd = currS // exclusive
val gapLen = gapEnd - gapStart
if (gapLen <= 0) {
// No gap → merge if length OK
if ((currE - prevS + 1) <= maxMergedLen) {
merged[merged.lastIndex] = prevS..currE
} else {
merged.add(seg)
}
continue
}
// C1: posterior wall ascent
if (x[currE] <= x[prevE]) { merged.add(seg); continue }
// C2: gap length
if (gapLen > maxGapLen) { merged.add(seg); continue }
// C3: gap median vs ref mean
val gap = DoubleArray(gapLen) { x[gapStart + it] }
val gapMedian = npMedian(gap)
val wPrev = min(refWin, prevE - prevS + 1)
val wCurr = min(refWin, currE - currS + 1)
var refSum = 0.0
for (k in 0 until wPrev) refSum += x[prevE - wPrev + 1 + k]
for (k in 0 until wCurr) refSum += x[currS + k]
val refMean = refSum / (wPrev + wCurr)
if (refMean == 0.0) { merged.add(seg); continue }
val relDiff = abs(gapMedian - refMean) / abs(refMean)
if (relDiff >= relTol) { merged.add(seg); continue }
// C4: merged length
if ((currE - prevS + 1) > maxMergedLen) { merged.add(seg); continue }
// merge
merged[merged.lastIndex] = prevS..currE
}
return merged
}
}
@@ -0,0 +1,147 @@
/*
* Generic Savitzky-Golay smoother — port of vesiscan_test/library/denoising.py:sg_smooth.
*
* 임의 (window, polyorder) 에 대해 Python 1:1 동작:
* 1) 내부 m..n-m: pinv(Vandermonde)[0] 커널로 컨볼루션
* 2) 왼쪽 edge 0..m-1: 첫 window 샘플에 polynomial fit → t=i-m 위치 평가
* 3) 오른쪽 edge n-m..n-1: 마지막 window 샘플에 polynomial fit → t=i-(n-m-1) 위치 평가
* 4) n < window: 전체 신호 단일 polynomial fit
*
* config_6ch.py: SG_WIN=7, SG_POLY=3 ← method_d 기준
* (V4.1 detector 는 별도 (5,2) 하드코딩 커널 사용 — walldetect/algo/Denoising.kt)
*
* 수치 검증:
* (7,3) 내부 커널 = [-2, 3, 6, 7, 6, 3, -2] / 21 (표준 SG 7-3 좌표)
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
object SgSmoothGeneric {
/**
* SG smooth signal `x` with given (window, polyorder).
* window 은 odd, polyorder < window 이어야 함.
*/
fun smooth(x: DoubleArray, window: Int, polyorder: Int): DoubleArray {
require(window % 2 == 1) { "window must be odd, got $window" }
require(polyorder < window) { "polyorder($polyorder) must be < window($window)" }
val n = x.size
val m = (window - 1) / 2
if (n < window) {
// 짧은 신호: 전체 구간 단일 polynomial fit
val xs = DoubleArray(n) { it - (n - 1) / 2.0 }
val coef = polyfit(xs, x, polyorder)
return DoubleArray(n) { i -> evalPoly(coef, xs[i]) }
}
val xsInner = DoubleArray(window) { it - m.toDouble() }
// 내부 커널 = pinv(A)[0, :] — 다항식 c0 (상수항) 의 LS 계수
val kernel = innerKernel(xsInner, polyorder)
val out = DoubleArray(n)
for (i in m until n - m) {
var s = 0.0
for (k in 0 until window) s += kernel[k] * x[i - m + k]
out[i] = s
}
// 왼쪽 edge: 첫 window 샘플 polynomial fit
val leftCoef = polyfit(xsInner, sliceArray(x, 0, window), polyorder)
for (i in 0 until m) {
val t = (i - m).toDouble() // block center(index m) 기준 상대좌표
out[i] = evalPoly(leftCoef, t)
}
// 오른쪽 edge: 마지막 window 샘플 polynomial fit
val rightCoef = polyfit(xsInner, sliceArray(x, n - window, n), polyorder)
for (i in n - m until n) {
val t = (i - (n - m - 1)).toDouble() // block center(index n-m-1) 기준 상대좌표
out[i] = evalPoly(rightCoef, t)
}
return out
}
private fun sliceArray(x: DoubleArray, from: Int, to: Int): DoubleArray =
DoubleArray(to - from) { x[from + it] }
private fun evalPoly(coef: DoubleArray, t: Double): Double {
var v = 0.0
var tk = 1.0
for (k in coef.indices) {
v += coef[k] * tk
tk *= t
}
return v
}
/** xs (length=window) 에 y (length=window) polynomial-fit → coefs [c0, c1, ..., c_p]. */
private fun polyfit(xs: DoubleArray, ys: DoubleArray, polyorder: Int): DoubleArray {
val p = polyorder + 1
val n = xs.size
val ata = Array(p) { DoubleArray(p) }
val aty = DoubleArray(p)
for (i in 0 until n) {
val powers = DoubleArray(p)
powers[0] = 1.0
for (k in 1 until p) powers[k] = powers[k - 1] * xs[i]
for (j in 0 until p) {
aty[j] += powers[j] * ys[i]
for (k in 0 until p) ata[j][k] += powers[j] * powers[k]
}
}
return solveLinearSystem(ata, aty)
}
/** 내부 SG 커널: pinv(A)[0, :] — c0 의 LS 계수. y → c0 = Σ kernel[i]·y[i]. */
private fun innerKernel(xs: DoubleArray, polyorder: Int): DoubleArray {
val p = polyorder + 1
val n = xs.size
val ata = Array(p) { DoubleArray(p) }
for (i in 0 until n) {
val powers = DoubleArray(p)
powers[0] = 1.0
for (k in 1 until p) powers[k] = powers[k - 1] * xs[i]
for (j in 0 until p) for (k in 0 until p) ata[j][k] += powers[j] * powers[k]
}
val e0 = DoubleArray(p)
e0[0] = 1.0
val nInvCol0 = solveLinearSystem(ata, e0) // N^-1 [:, 0] = 대칭으로 row 0 동등
val kernel = DoubleArray(n)
for (i in 0 until n) {
var xik = 1.0
for (k in 0 until p) {
kernel[i] += nInvCol0[k] * xik
xik *= xs[i]
}
}
return kernel
}
/** Gaussian elimination with partial pivoting — 소형 PSD 정상 매트릭스용. */
private fun solveLinearSystem(a: Array<DoubleArray>, b: DoubleArray): DoubleArray {
val n = b.size
val m = Array(n) { DoubleArray(n + 1) }
for (i in 0 until n) {
for (j in 0 until n) m[i][j] = a[i][j]
m[i][n] = b[i]
}
for (i in 0 until n) {
var pivot = i
for (k in i + 1 until n) if (abs(m[k][i]) > abs(m[pivot][i])) pivot = k
if (pivot != i) { val t = m[i]; m[i] = m[pivot]; m[pivot] = t }
val piv = m[i][i]
require(abs(piv) >= 1e-12) { "singular matrix at row $i" }
for (k in i + 1 until n) {
val f = m[k][i] / piv
for (j in i..n) m[k][j] -= f * m[i][j]
}
}
val x = DoubleArray(n)
for (i in n - 1 downTo 0) {
var s = m[i][n]
for (j in i + 1 until n) s -= m[i][j] * x[j]
x[i] = s / m[i][i]
}
return x
}
}
@@ -0,0 +1,77 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/single_channel_bv.js (1:1).
*
* Chord-based bladder volume from a single channel.
* Sphere assumption: chord ≡ 2R, so BV = (4/3)π(c/2)³.
* §11 Mode B fallback for the small-bladder regime.
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdConfig
import kotlin.math.cbrt
object SingleChannelBv {
enum class Reduction { MEDIAN, MEAN }
data class PerChannel(
val ch: Int,
val ant: Int?,
val post: Int?,
val chord: Double?, // mm
val r: Double?, // mm (chord / 2)
val bv: Double? // mL
)
data class AggResult(
val bv: Double?,
val r: Double?,
val n: Int,
val used: List<PerChannel>
)
fun chordBVPerChannel(
detResults: List<Geometry.Detection>,
dpsMm: Double = WdConfig.DPS_DEFAULT
): List<PerChannel> {
return detResults.mapIndexed { ch, r ->
if (r.ant == null || r.post == null) {
PerChannel(ch = ch, ant = null, post = null, chord = null, r = null, bv = null)
} else {
val chord = (r.post - r.ant) * dpsMm
val rMm = chord / 2.0
val bv = (4.0 / 3.0) * Math.PI * rMm * rMm * rMm / 1000.0
PerChannel(ch = ch, ant = r.ant, post = r.post, chord = chord, r = rMm, bv = bv)
}
}
}
/**
* Aggregate single-channel BVs from trusted channels.
* @param trusted ch → trust flag (already filtered by s_contrast)
* @param reduction MEDIAN (default — JS upper-median floor(N/2)) or MEAN
*/
fun aggregate(
perChannelBV: List<PerChannel>,
trusted: Map<Int, Boolean>,
reduction: Reduction = Reduction.MEDIAN
): AggResult {
val used = perChannelBV.filter { it.bv != null && trusted[it.ch] == true }
val vals = used.mapNotNull { it.bv }
if (vals.isEmpty()) return AggResult(bv = null, r = null, n = 0, used = emptyList())
val bv = if (reduction == Reduction.MEAN) {
vals.average()
} else {
// JS upper-median: sorted[floor(N/2)]
val sorted = vals.sorted()
sorted[sorted.size / 2]
}
val rMm = cbrt(3.0 * bv * 1000.0 / (4.0 * Math.PI))
return AggResult(bv = bv, r = rMm, n = vals.size, used = used)
}
}
@@ -0,0 +1,124 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/span_utils.js (1:1).
* Contiguous boolean spans + close-span merge with peak guard.
*/
package com.medithings.vesiscan.walldetect.algo
object SpanUtils {
/** A boolean-mask span [start, end] inclusive. Equivalent to JS [s, e] tuple. */
data class Span(val start: Int, val end: Int)
fun contiguousTrueSpans(mask: BooleanArray): List<Span> {
val spans = mutableListOf<Span>()
var start = -1
for (i in mask.indices) {
if (mask[i] && start == -1) {
start = i
} else if (!mask[i] && start != -1) {
spans += Span(start, i - 1)
start = -1
}
}
if (start != -1) spans += Span(start, mask.size - 1)
return spans
}
/**
* Stage 1 merge — width-bounded amplitude-aware merge. Closes gaps
* of width ≤ maxGap when the gap-peak amplitude is below gapPeakThr
* (typically T + 50 ADC). Compatible with V2 / Otsu / scalar paths.
*
* Note: gap-peak should be sampled from the ORIGINAL `sg` (not the
* median-filtered envelope) because the running median can clip a
* real wall peak to its plateau-median value, which on borderline
* cases drops below T + GAP_PEAK_MARGIN and would erroneously merge
* across the wall.
*/
fun mergeCloseSpans(
spans: List<Span>,
maxGap: Int,
sg: DoubleArray? = null,
gapPeakThr: Double? = null
): List<Span> {
if (spans.isEmpty()) return emptyList()
val ordered = spans.sortedBy { it.start }
val merged = mutableListOf(ordered[0])
val usePeak = sg != null && gapPeakThr != null
for (k in 1 until ordered.size) {
val (s, e) = ordered[k]
val last = merged.last()
val pe = last.end
val gap = s - pe - 1
if (gap <= maxGap) {
var skip = false
if (usePeak && gap >= 1) {
var mx = Double.NEGATIVE_INFINITY
for (i in (pe + 1) until s) if (sg!![i] > mx) mx = sg[i]
if (mx > gapPeakThr!!) skip = true
}
if (!skip) {
merged[merged.lastIndex] = Span(last.start, maxOf(pe, e))
continue
}
}
merged += Span(s, e)
}
return merged
}
/**
* Stage 2 bimodal merge (V4.1 only) — Otsu 1979 split between
* lumen-baseline and wall-echo classes, with an anatomical
* postMaxIdx ceiling. No width cap (handles wide intra-lumen
* speckle clusters that exceed the median's ⌊W/2⌋ absorption
* width), but the merged span end must remain inside the
* plausible-wall depth range so that real post-wall tails are not
* absorbed.
*
* Caller should gate this stage on `spans.size >= 3` (a normal
* capture leaves stage 1 with exactly 2 spans — lumen + tail —
* and needs no further merging).
*
* §3.7 Dual-threshold guard (added 2026-05): when `gapPeakHi` is
* provided (typically `T_cfar + GAP_PEAK_MARGIN_STAGE2`), the gap
* peak must fall BELOW BOTH `otsuThr` AND `gapPeakHi`. This blocks
* the wall+container double-peak failure mode (Japan-standard
* 150 mL CH0/CH5: a strong reflector beyond the bladder pulls
* Otsu's threshold up, causing the legitimate intermediate wall
* echo to be mis-classified as speckle). The CFAR-derived ceiling
* is calibrated to lumen-noise statistics; the AND combination
* provides cross-validation across two orthogonal histograms.
*/
fun mergeBimodal(
spans: List<Span>,
sg: DoubleArray,
otsuThr: Double,
postMaxIdx: Int,
gapPeakHi: Double? = null
): List<Span> {
if (spans.isEmpty()) return emptyList()
val ordered = spans.sortedBy { it.start }
val merged = mutableListOf(ordered[0])
for (k in 1 until ordered.size) {
val (s, e) = ordered[k]
val last = merged.last()
val pe = last.end
var mx = Double.NEGATIVE_INFINITY
for (i in (pe + 1) until s) if (sg[i] > mx) mx = sg[i]
val newEnd = maxOf(pe, e)
val passOtsu = mx < otsuThr
val passGuard = (gapPeakHi == null) || (mx < gapPeakHi)
if (passOtsu && passGuard && newEnd <= postMaxIdx) {
merged[merged.lastIndex] = Span(last.start, newEnd)
} else {
merged += Span(s, e)
}
}
return merged
}
}
@@ -0,0 +1,159 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/sphere_fit_2step.js (1:1).
*
* Kasa → LM, the canonical CharlesKWON V4 sphere fit.
* Auto-selects circle (xz / yz / xy plane) when all points share one axis,
* else 3D sphere.
*
* 'auto' detection (1e-6 epsilon):
* constY → circle (drop y, xz plane) — V2 30° SI-only standard case
* constX → circle (drop x, yz plane) — central-column-only post-gate
* constZ → circle (drop z, xy plane)
* else → sphere (3D)
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
import kotlin.math.sqrt
object SphereFit2Step {
enum class Mode { AUTO, CIRCLE, SPHERE }
data class ResidualStats(
val residuals: DoubleArray, // (|P_i − C| − R)
val mean: Double,
val std: Double,
val max: Double // max(|residual|)
)
/**
* Fit result. `mode` is the actually-chosen mode (auto resolves to circle/sphere).
* `dropAxis` set only for circle: 0=x, 1=y, 2=z.
* `axes` is the index pair of the spanning plane (only set for circle).
* `kasaCenter`/`kasaR` is the initial guess; `lmCenter`/`lmR` is the refined fit.
* `residualStats` is computed against the LM solution.
*/
data class FitResult(
val mode: String, // "circle" | "sphere"
val dropAxis: Int?, // 0/1/2 for circle, null for sphere
val axes: IntArray?, // [a, b] spanning plane (circle only), null for sphere
val kasaCenter: DoubleArray,
val kasaR: Double,
val lmCenter: DoubleArray, // 2 entries for circle, 3 for sphere
val lmR: Double,
val lmIter: Int,
val residuals: DoubleArray,
val residualMean: Double,
val residualStd: Double,
val residualMax: Double
)
/** Distance from each point to centre. */
fun residuals(p: List<DoubleArray>, center: DoubleArray): DoubleArray {
val k = center.size
return DoubleArray(p.size) { i ->
var s = 0.0
for (d in 0 until k) {
val dd = p[i][d] - center[d]
s += dd * dd
}
sqrt(s)
}
}
fun residualStats(p: List<DoubleArray>, center: DoubleArray, r: Double): ResidualStats {
val raw = residuals(p, center)
val res = DoubleArray(raw.size) { raw[it] - r }
var sum = 0.0
for (v in res) sum += v
val mean = sum / res.size
var sq = 0.0
for (v in res) { val d = v - mean; sq += d * d }
var maxAbs = 0.0
for (v in res) { val a = abs(v); if (a > maxAbs) maxAbs = a }
return ResidualStats(
residuals = res,
mean = mean,
std = sqrt(sq / res.size),
max = maxAbs
)
}
/** Project P (N×3) onto the 2D plane spanned by axes (a, b). */
private fun project(p: List<DoubleArray>, a: Int, b: Int): List<DoubleArray> =
p.map { doubleArrayOf(it[a], it[b]) }
fun fit2Step(p: List<DoubleArray>, mode: Mode = Mode.AUTO): FitResult? {
if (p.size < 3) return null
val is3D = p.all { it.size == 3 }
var resolvedMode = mode
var dropAxis: Int? = null
if (mode == Mode.AUTO && is3D) {
val eps = 1e-6
val constX = p.all { abs(it[0] - p[0][0]) < eps }
val constY = p.all { abs(it[1] - p[0][1]) < eps }
val constZ = p.all { abs(it[2] - p[0][2]) < eps }
when {
constY -> { resolvedMode = Mode.CIRCLE; dropAxis = 1 } // historical V2 30°
constX -> { resolvedMode = Mode.CIRCLE; dropAxis = 0 } // central column only
constZ -> { resolvedMode = Mode.CIRCLE; dropAxis = 2 }
else -> { resolvedMode = Mode.SPHERE }
}
} else if (mode == Mode.CIRCLE && is3D && dropAxis == null) {
// Caller forced circle without indicating axis → historical default (drop y, fit xz)
dropAxis = 1
}
return if (resolvedMode == Mode.CIRCLE) {
val axes = when (dropAxis) {
0 -> intArrayOf(1, 2)
1 -> intArrayOf(0, 2)
2 -> intArrayOf(0, 1)
else -> intArrayOf(0, 1) // 2D input direct; axes meaningless
}
val pPlane = if (is3D) project(p, axes[0], axes[1]) else p
val k0 = SphereKasa.kasaCircle(pPlane) ?: return null
val fin = SphereLm.lmCircle(pPlane, k0.center, k0.r)
val stats = residualStats(pPlane, fin.center, fin.r)
FitResult(
mode = "circle",
dropAxis = dropAxis,
axes = axes,
kasaCenter = k0.center,
kasaR = k0.r,
lmCenter = fin.center,
lmR = fin.r,
lmIter = fin.iter,
residuals = stats.residuals,
residualMean = stats.mean,
residualStd = stats.std,
residualMax = stats.max
)
} else {
val k0 = SphereKasa.kasaSphere(p) ?: return null
val fin = SphereLm.lmSphere(p, k0.center, k0.r)
val stats = residualStats(p, fin.center, fin.r)
FitResult(
mode = "sphere",
dropAxis = null,
axes = null,
kasaCenter = k0.center,
kasaR = k0.r,
lmCenter = fin.center,
lmR = fin.r,
lmIter = fin.iter,
residuals = stats.residuals,
residualMean = stats.mean,
residualStd = stats.std,
residualMax = stats.max
)
}
}
}
@@ -0,0 +1,60 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/sphere_kasa.js (1:1).
*
* Kasa (1976) algebraic least-squares circle/sphere fit. Linearises
* |p|² = 2 c·p + (R² − |c|²)
* with unknowns (c, R²−|c|²). Used as the Mode A initial guess before LM.
*
* Reference:
* Kasa I., "A circle fitting procedure and its error analysis",
* IEEE Trans Instrum Meas IM-25(1):8–14, 1976.
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdLinalg
import kotlin.math.max
import kotlin.math.sqrt
object SphereKasa {
data class CircleFit(val center: DoubleArray, val r: Double)
data class SphereFit(val center: DoubleArray, val r: Double)
/** P : list of (x, y) pairs. Returns null if n < 3 or singular. */
fun kasaCircle(p: List<DoubleArray>): CircleFit? {
val n = p.size
if (n < 3) return null
val a = Array(n) { DoubleArray(3) }
val b = DoubleArray(n)
for (i in 0 until n) {
val x = p[i][0]; val y = p[i][1]
a[i][0] = 2 * x; a[i][1] = 2 * y; a[i][2] = 1.0
b[i] = x * x + y * y
}
val sol = WdLinalg.solve(WdLinalg.ata(a), WdLinalg.atb(a, b)) ?: return null
val cx = sol[0]; val cy = sol[1]; val c = sol[2]
val r2 = c + cx * cx + cy * cy
return CircleFit(center = doubleArrayOf(cx, cy), r = sqrt(max(r2, 0.0)))
}
/** P : list of (x, y, z) triples. Returns null if n < 4 or singular. */
fun kasaSphere(p: List<DoubleArray>): SphereFit? {
val n = p.size
if (n < 4) return null
val a = Array(n) { DoubleArray(4) }
val b = DoubleArray(n)
for (i in 0 until n) {
val x = p[i][0]; val y = p[i][1]; val z = p[i][2]
a[i][0] = 2 * x; a[i][1] = 2 * y; a[i][2] = 2 * z; a[i][3] = 1.0
b[i] = x * x + y * y + z * z
}
val sol = WdLinalg.solve(WdLinalg.ata(a), WdLinalg.atb(a, b)) ?: return null
val cx = sol[0]; val cy = sol[1]; val cz = sol[2]; val c = sol[3]
val r2 = c + cx * cx + cy * cy + cz * cz
return SphereFit(center = doubleArrayOf(cx, cy, cz), r = sqrt(max(r2, 0.0)))
}
}
@@ -0,0 +1,95 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/sphere_lm.js (1:1).
*
* Levenberg-Marquardt geometric LS for circle / sphere.
* Minimise Σ (|P_i − C| − R)²
* Jacobian: ∂r_i/∂c = −(P_i − C) / |P_i − C|; ∂r_i/∂R = −1
*
* Reference:
* Levenberg K. (1944) Quart. Appl. Math. 2:164–168.
* Marquardt D. (1963) SIAM J. Appl. Math. 11:431–441.
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdLinalg
import kotlin.math.sqrt
object SphereLm {
data class CircleResult(val center: DoubleArray, val r: Double, val iter: Int)
data class SphereResult(val center: DoubleArray, val r: Double, val iter: Int)
fun lmCircle(
p: List<DoubleArray>,
c0: DoubleArray,
r0: Double,
maxIter: Int = 50,
tol: Double = 1e-8,
lam: Double = 1e-3
): CircleResult {
var cx = c0[0]; var cy = c0[1]; var r = r0
var iter = 0
while (iter < maxIter) {
val n = p.size
val j = Array(n) { DoubleArray(3) }
val res = DoubleArray(n)
for (i in 0 until n) {
val dx = p[i][0] - cx
val dy = p[i][1] - cy
val d = sqrt(dx * dx + dy * dy).let { if (it == 0.0) 1e-12 else it }
j[i][0] = -dx / d
j[i][1] = -dy / d
j[i][2] = -1.0
res[i] = d - r
}
val h = WdLinalg.ata(j)
h[0][0] += lam; h[1][1] += lam; h[2][2] += lam
val g = WdLinalg.atb(j, res)
val delta = WdLinalg.solve(h, doubleArrayOf(-g[0], -g[1], -g[2])) ?: break
cx += delta[0]; cy += delta[1]; r += delta[2]
iter++
if (WdLinalg.vecNorm(delta) < tol) break
}
return CircleResult(center = doubleArrayOf(cx, cy), r = r, iter = iter + 1)
}
fun lmSphere(
p: List<DoubleArray>,
c0: DoubleArray,
r0: Double,
maxIter: Int = 50,
tol: Double = 1e-8,
lam: Double = 1e-3
): SphereResult {
var cx = c0[0]; var cy = c0[1]; var cz = c0[2]; var r = r0
var iter = 0
while (iter < maxIter) {
val n = p.size
val j = Array(n) { DoubleArray(4) }
val res = DoubleArray(n)
for (i in 0 until n) {
val dx = p[i][0] - cx
val dy = p[i][1] - cy
val dz = p[i][2] - cz
val d = sqrt(dx * dx + dy * dy + dz * dz).let { if (it == 0.0) 1e-12 else it }
j[i][0] = -dx / d
j[i][1] = -dy / d
j[i][2] = -dz / d
j[i][3] = -1.0
res[i] = d - r
}
val h = WdLinalg.ata(j)
for (k in 0..3) h[k][k] += lam
val g = WdLinalg.atb(j, res)
val delta = WdLinalg.solve(h, doubleArrayOf(-g[0], -g[1], -g[2], -g[3])) ?: break
cx += delta[0]; cy += delta[1]; cz += delta[2]; r += delta[3]
iter++
if (WdLinalg.vecNorm(delta) < tol) break
}
return SphereResult(center = doubleArrayOf(cx, cy, cz), r = r, iter = iter + 1)
}
}
@@ -0,0 +1,88 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/sta_lta.js (1:1).
*
* Short-Term Average / Long-Term Average impulse detector — Allen 1978.
*
* Reference (foundational)
* Allen RV. "Automatic earthquake recognition and timing from single
* traces." Bull Seismol Soc Am 68(5):1521-1532, 1978.
*
* Trnkoczy A. "Understanding and parameter setting of STA/LTA trigger
* algorithm." in IASPEI New Manual of Seismological Observatory
* Practice (NMSOP-2) §8.1, 2012. doi:10.2312/GFZ.NMSOP-2_IS_8.1
*
* Withers M, Aster R, Young C, et al. "A comparison of select trigger
* algorithms for automated global seismic phase and event detection."
* Bull Seismol Soc Am 88(1):95-106, 1998.
*
* Used in V4.1 ONLY by BModeScore for the impulse-purity subscore u_stl —
* discriminates the sharp wall+floor merged echo of phantom-on-rigid-floor
* captures from smooth reverberation bumps.
*
* STA[i] = (1/Nsta)·Σ_{k=i-Nsta+1..i} r²[k]
* LTA[i] = (1/Nlta)·Σ_{k=i-Nlta+1..i} r²[k]
* R[i] = STA[i] / LTA[i]
*
* Defaults (Trnkoczy 2012 §8.1.2): Nsta=3, Nlta=30 — tuned for short-
* duration impulse (1-3 samples) in stationary background noise.
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.max
import kotlin.math.min
object StaLta {
const val DEFAULT_NSTA = 3
const val DEFAULT_NLTA = 30
data class Result(val ratio: DoubleArray, val sta: DoubleArray, val lta: DoubleArray)
/**
* Full per-sample STA/LTA ratio of envelope².
* The first (Nlta-1) samples have ratio set to 1.0 (no LTA history yet)
* to suppress spurious early triggers.
*/
fun compute(envelope: DoubleArray, nSta: Int = DEFAULT_NSTA, nLta: Int = DEFAULT_NLTA): Result {
val n = envelope.size
val sq = DoubleArray(n) { envelope[it] * envelope[it] }
val cum = DoubleArray(n + 1)
for (i in 0 until n) cum[i + 1] = cum[i] + sq[i]
val sta = DoubleArray(n)
val lta = DoubleArray(n)
val ratio = DoubleArray(n)
for (i in 0 until n) {
val sLo = max(0, i - nSta + 1)
val sHi = i + 1
sta[i] = (cum[sHi] - cum[sLo]) / (sHi - sLo)
val lLo = max(0, i - nLta + 1)
val lHi = i + 1
lta[i] = (cum[lHi] - cum[lLo]) / (lHi - lLo)
ratio[i] = if (i < nLta - 1) 1.0
else if (lta[i] > 0.0) sta[i] / lta[i] else 0.0
}
return Result(ratio, sta, lta)
}
/**
* Convenience: maximum STA/LTA ratio in a small window around `idx`
* ([idx-2, idx+5]). Used to estimate impulse purity at a known wall
* position (V4.1 BModeScore.u_stl subscore).
*/
fun peakRatio(envelope: DoubleArray, idx: Int?, nSta: Int = DEFAULT_NSTA, nLta: Int = DEFAULT_NLTA): Double {
if (idx == null) return 0.0
val n = envelope.size
if (n == 0) return 0.0
val lo = max(0, idx - 2)
val hi = min(n - 1, idx + 5)
val r = compute(envelope, nSta, nLta).ratio
var m = 0.0
for (i in lo..hi) if (r[i] > m) m = r[i]
return m
}
}
@@ -0,0 +1,81 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/subsample_refine.js (1:1).
* Sub-sample wall position refinement (Cespedes 1995, parabolic 3-point fit).
*
* Method:
* y(x) = y0 + ½ y'' (x − x0)²
* 3-point fit (idx-1, idx, idx+1) gives parabola vertex at
* δ = ½ · (y[idx-1] − y[idx+1]) / (y[idx-1] − 2·y[idx] + y[idx+1])
* |δ| ≤ 0.5 when idx is a strict local extremum (else fall back to integer idx).
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
object SubsampleRefine {
enum class Kind { AUTO, PEAK, VALLEY }
/**
* @return refined fractional sample index, or null if `idx` is null,
* or the integer `idx` if refinement is not applicable.
*/
fun refineParabolic(envelope: DoubleArray?, idx: Int?, kind: Kind = Kind.AUTO): Double? {
if (envelope == null || idx == null) return idx?.toDouble()
val n = envelope.size
if (idx < 1 || idx > n - 2) return idx.toDouble()
val ym1 = envelope[idx - 1]
val y0 = envelope[idx]
val yp1 = envelope[idx + 1]
val denom = ym1 - 2 * y0 + yp1
if (denom == 0.0 || !denom.isFinite()) return idx.toDouble()
if (kind == Kind.PEAK && denom > 0) return idx.toDouble()
if (kind == Kind.VALLEY && denom < 0) return idx.toDouble()
val delta = 0.5 * (ym1 - yp1) / denom
if (!delta.isFinite() || abs(delta) > 1.0) return idx.toDouble()
return idx + delta
}
/**
* Linear interpolation of an envelope crossing T near `idx`.
* @param dir "up" — sign change neg → non-neg (ant-like)
* "down" — sign change non-neg → neg (post-like)
* null — accept either direction
*/
fun refineLinearCrossing(
envelope: DoubleArray?,
idx: Int?,
T: Double,
dir: String?
): Double? {
if (envelope == null || idx == null) return idx?.toDouble()
val n = envelope.size
if (idx < 1 || idx > n - 1) return idx.toDouble()
val samples = listOf(idx - 1, idx, idx + 1).filter { it in 0 until n }
for (s in 0 until samples.size - 1) {
val a = samples[s]
val b = samples[s + 1]
val va = envelope[a] - T
val vb = envelope[b] - T
val isUp = (va < 0 && vb >= 0)
val isDown = (va >= 0 && vb < 0)
if ((dir == "up" && isUp) || (dir == "down" && isDown) ||
(dir == null && (isUp || isDown))
) {
val denom = vb - va
if (denom == 0.0) return a.toDouble()
val frac = -va / denom // ∈ [0, 1]
return a + frac
}
}
return idx.toDouble()
}
}
@@ -0,0 +1,72 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/threshold_oscfar.js (1:1).
*
* Order-Statistic CFAR (Rohling 1983) — adaptive threshold in V4.1.
* T_i = scale · Q_rank(window_i \ {i}) [per-sample]
* T = median(T_i) [scalar]
*
* Reference:
* Rohling H., "Radar CFAR Thresholding in Clutter and Multiple Target
* Situations", IEEE Trans Aerosp Electron Syst AES-19(4):608–621, 1983.
*
* Defaults (CONFIG):
* OSCFAR_WIN = 15 (centred window length, half = 7 each side)
* OSCFAR_RANK = 0.4 (40th percentile in window-without-CUT)
* OSCFAR_SCALE = 1.05 (5% safety margin above the rank percentile)
*/
package com.medithings.vesiscan.walldetect.algo
import com.medithings.vesiscan.walldetect.core.WdConfig
import com.medithings.vesiscan.walldetect.core.WdNumeric
import kotlin.math.max
import kotlin.math.min
object ThresholdOsCfar {
/**
* Per-sample OS-CFAR threshold trace.
* For each sample i: collect window samples (i ± half) excluding i itself,
* compute the rank-th quantile, scale by `scale`. Edge windows are truncated.
* If the local set is empty (n=1 input), threshold = sg[i].
*/
fun perSample(
sg: DoubleArray,
win: Int = WdConfig.OSCFAR_WIN,
rank: Double = WdConfig.OSCFAR_RANK,
scale: Double = WdConfig.OSCFAR_SCALE
): DoubleArray {
val n = sg.size
val half = win / 2
val thr = DoubleArray(n)
for (i in 0 until n) {
val lo = max(0, i - half)
val hi = min(n, i + half + 1)
val localSize = (i - lo) + (hi - (i + 1))
if (localSize <= 0) {
thr[i] = sg[i]
continue
}
val local = DoubleArray(localSize)
var p = 0
for (k in lo until i) { local[p++] = sg[k] }
for (k in (i + 1) until hi) { local[p++] = sg[k] }
thr[i] = WdNumeric.quantile(local, rank) * scale
}
return thr
}
/**
* Scalar OS-CFAR threshold = median of per-sample trace.
* This is what V4.1 lumen-first detector uses as the lumen mask cut value.
*/
fun scalar(
sg: DoubleArray,
win: Int = WdConfig.OSCFAR_WIN,
rank: Double = WdConfig.OSCFAR_RANK,
scale: Double = WdConfig.OSCFAR_SCALE
): Double = WdNumeric.median(perSample(sg, win, rank, scale))
}
@@ -0,0 +1,132 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/wall_select.js (1:1).
* V2 prominence-based wall selection with other_edge inner clamp.
* Mirror of py2/low_echo_detection_method_b._select_wall_by_prominence.
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
import kotlin.math.max
import kotlin.math.min
object WallSelect {
const val PEAK_SEARCH_WIN = 20
/**
* Default candidate cap (legacy). Used by V2 and any caller that
* doesn't pass an explicit maxCandidates. V4.1 uses the side-aware
* defaults below.
*/
const val MAX_PEAK_CANDIDATES = 3
/**
* V4.1 ANT side — keep "closest 3" semantics. Anterior wall must be
* the LAST prominent peak just before the lumen begins. Letting
* prominence-only pick freely would elect a far-away transducer
* ring-down peak whose vertical depth then falls below the
* anatomical-gate floor.
*/
const val MAX_PEAK_CANDIDATES_ANT = 3
/**
* V4.1 POST side — large cap. On phantom-on-rigid-floor captures
* the bladder posterior wall + container floor merge into a single
* dominant peak that can be 5-15 samples FARTHER from the lumen
* edge than smaller intra-tissue ripples. Top-3-closest excludes
* it. Lifting the cap to 64 lets prominence — exactly the right
* discriminator — actually decide. peakMin still removes noise.
*/
const val MAX_PEAK_CANDIDATES_POST = 64
enum class Side { ANT, POST }
/**
* @param sg envelope (smoothed)
* @param edge span endpoint (s for ant, e for post)
* @param searchWin search window size around edge
* @param peakMin minimum amplitude to consider a peak candidate
* @param side ANT → prominence vs right valley (urine direction)
* POST → prominence vs left valley (urine direction)
* @param otherEdge opposite span endpoint — inner search must not cross it.
* @param maxCandidates limit per side
* @return wall sample index, or null if no candidate
*/
fun selectWallByProminence(
sg: DoubleArray,
edge: Int,
searchWin: Int = PEAK_SEARCH_WIN,
peakMin: Double,
side: Side,
otherEdge: Int? = null,
maxCandidates: Int = MAX_PEAK_CANDIDATES,
edgeDistDecay: Double = 0.12, // py2 method_b EDGE_DIST_DECAY
valleyStopRise: Double = 50.0, // py2 method_b VALLEY_STOP_RISE
): Int? {
val n = sg.size
var leftLo: Int
var rightHi: Int
if (side == Side.ANT) {
leftLo = max(0, edge - searchWin)
rightHi = min(n - 1, edge + searchWin)
if (otherEdge != null) rightHi = min(rightHi, otherEdge)
} else { // POST
leftLo = max(0, edge - searchWin)
if (otherEdge != null) leftLo = max(leftLo, otherEdge)
rightHi = min(n - 1, edge + searchWin)
}
// py2 method_b — left slice extended to edge+1 so that an edge−1 peak
// is detectable (find_peaks_1d excludes the rightmost sample).
var leftCand: List<Int> = emptyList()
if (edge > leftLo) {
val leftEnd = min(edge + 1, n)
val seg = DoubleArray(leftEnd - leftLo) { sg[leftLo + it] }
val pks = PeakDetection.findPeaks1D(seg)
leftCand = pks.map { leftLo + it }
.filter { it < edge && sg[it] >= peakMin }
.sortedBy { abs(it - edge) }
.take(maxCandidates)
}
// py2 method_b — right slice starts at edge−1 so that an edge peak is detectable.
var rightCand: List<Int> = emptyList()
if (rightHi >= edge) {
val rightStart = max(edge - 1, 0)
val seg = DoubleArray(rightHi + 1 - rightStart) { sg[rightStart + it] }
val pks = PeakDetection.findPeaks1D(seg)
rightCand = pks.map { rightStart + it }
.filter { it >= edge && sg[it] >= peakMin }
.sortedBy { abs(it - edge) }
.take(maxCandidates)
}
val candidates = leftCand + rightCand
if (candidates.isEmpty()) return null
// py2 method_b prominence + edge-distance decay:
// score = prominence / (1 + EDGE_DIST_DECAY * |p - edge|)
// Valley walk uses VALLEY_STOP_RISE instead of legacy 10.
var best: Int? = null
var bestScore = Double.NEGATIVE_INFINITY
for (p in candidates) {
val valley = if (side == Side.ANT)
PeakDetection.rightValley(sg, p, breakRise = valleyStopRise)
else
PeakDetection.leftValley(sg, p, breakRise = valleyStopRise)
val prom = sg[p] - valley
val dist = abs(p - edge).toDouble()
val score = prom / (1.0 + edgeDistDecay * dist)
if (score > bestScore) {
bestScore = score
best = p
}
}
return best
}
}
@@ -0,0 +1,309 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/algo/wavelet_denoise.js (1:1).
*
* Sun 2024 wavelet energy-ratio adaptive denoising for A-mode envelopes.
* Reference: Sun J et al., "A-Mode Ultrasound Bladder Volume Estimation
* Algorithm Based on Wavelet Energy Ratio Adaptive Denoising,"
* Sensors 24(6):1984, 2024. doi:10.3390/s24061984
*
* Pipeline:
* 1. Pad signal to multiple of 2^L (periodic extension at end).
* 2. L-level Daubechies-4 DWT decomposition.
* 3. Per-level noise sigma_j = MAD(d_j) × 1.4826 (Donoho 1995).
* 4. Per-level energy E_j = ||d_j||² / N_j; ratio R_j = E_j / max_k(E_k).
* 5. Adaptive soft threshold: λ_j = sigma_j · √(2·ln N_j) · ((1 − R_j) + ε).
* 6. Inverse DWT, truncate to original length.
*
* Convention: details[0] = level 1 (highest freq), details[L-1] = level L
* (lowest detail / nearest to approx). approx is stored separately.
*/
package com.medithings.vesiscan.walldetect.algo
import kotlin.math.abs
import kotlin.math.ceil
import kotlin.math.floor
import kotlin.math.ln
import kotlin.math.max
import kotlin.math.min
import kotlin.math.sqrt
object WaveletDenoise {
// ── DB4 (Daubechies-4) filter coefficients — PyWavelets canonical ──
val DB4_DEC_LO = doubleArrayOf(
-0.010597401784997278, 0.032883011666982945, 0.030841381835986965, -0.18703481171888114,
-0.027983769416983849, 0.63088076792959036, 0.71484657055254153, 0.23037781330885523
)
val DB4_DEC_HI = doubleArrayOf(
-0.23037781330885523, 0.71484657055254153, -0.63088076792959036, -0.027983769416983849,
0.18703481171888114, 0.030841381835986965, -0.032883011666982945, -0.010597401784997278
)
val DB4_REC_LO = doubleArrayOf(
0.23037781330885523, 0.71484657055254153, 0.63088076792959036, -0.027983769416983849,
-0.18703481171888114, 0.030841381835986965, 0.032883011666982945, -0.010597401784997278
)
val DB4_REC_HI = doubleArrayOf(
-0.010597401784997278, -0.032883011666982945, 0.030841381835986965, 0.18703481171888114,
-0.027983769416983849, -0.63088076792959036, 0.71484657055254153, -0.23037781330885523
)
/** Single-level DWT result: a = approximation, d = detail. */
data class Dwt1Result(val a: DoubleArray, val d: DoubleArray)
/**
* Multi-level Mallat decomposition.
* details[0] = level 1 (highest freq), details[levels-1] = level L (deepest).
*/
data class Decomposition(
val approx: DoubleArray,
val details: List<DoubleArray>,
val lengths: IntArray,
val paddedLength: Int
)
/** Per-level diagnostics (energies, ratios, sigma estimates). */
data class Diagnose(
val levels: Int,
val sigma: DoubleArray,
val energy: DoubleArray,
val ratio: DoubleArray,
val detailLengths: IntArray
)
// ─────────────────────────────────────────────────────────
// Single-level DWT / IDWT (periodic boundary)
// ─────────────────────────────────────────────────────────
/** Standard convolution DWT: a[k] = Σ h_lo[m] · x[(2k − m + N) mod N]. */
fun dwt1(x: DoubleArray, hLo: DoubleArray, hHi: DoubleArray): Dwt1Result {
val n = x.size
val m = n shr 1
val l = hLo.size
val a = DoubleArray(m)
val d = DoubleArray(m)
for (k in 0 until m) {
var sa = 0.0
var sd = 0.0
for (nn in 0 until l) {
val i = (((2 * k - nn) % n) + n) % n
sa += hLo[nn] * x[i]
sd += hHi[nn] * x[i]
}
a[k] = sa
d[k] = sd
}
return Dwt1Result(a, d)
}
/**
* Synthesis with (L-1) sample shift to align idwt(dwt(x)) = x.
* out[n] = Σ_k g_lo[((n + L − 1) − 2k) mod N] · a[k]
* + Σ_k g_hi[((n + L − 1) − 2k) mod N] · d[k]
*/
fun idwt1(a: DoubleArray, d: DoubleArray, gLo: DoubleArray, gHi: DoubleArray): DoubleArray {
val mLen = a.size
val n = 2 * mLen
val l = gLo.size
val out = DoubleArray(n)
val shift = l - 1
for (nn in 0 until n) {
var s = 0.0
for (k in 0 until mLen) {
val mIdx = (((nn + shift) - 2 * k) % n + n) % n
if (mIdx < l) s += gLo[mIdx] * a[k] + gHi[mIdx] * d[k]
}
out[nn] = s
}
return out
}
// ─────────────────────────────────────────────────────────
// Multi-level Mallat pyramidal DWT/IDWT
// ─────────────────────────────────────────────────────────
fun dwt(signal: DoubleArray, levels: Int = 3): Decomposition {
val n0 = signal.size
val pow = 1 shl levels
val target = (ceil(n0.toDouble() / pow) * pow).toInt()
val x: DoubleArray = if (target != n0) {
DoubleArray(target) { signal[it % n0] }
} else {
signal.copyOf()
}
val details = mutableListOf<DoubleArray>()
val lengths = mutableListOf(x.size)
var approx = x
for (j in 0 until levels) {
val out = dwt1(approx, DB4_DEC_LO, DB4_DEC_HI)
details += out.d
approx = out.a
lengths += approx.size
}
return Decomposition(approx, details.toList(), lengths.toIntArray(), x.size)
}
fun idwt(decomp: Decomposition): DoubleArray {
var approx = decomp.approx
for (j in decomp.details.size - 1 downTo 0) {
approx = idwt1(approx, decomp.details[j], DB4_REC_LO, DB4_REC_HI)
}
return approx
}
// ─────────────────────────────────────────────────────────
// MAD (median absolute deviation) — JS uses simple median, NOT linear-interp
// ─────────────────────────────────────────────────────────
/** Mirror of JS `mad()`: median via floor(N/2), no linear interpolation. */
fun mad(arr: DoubleArray): Double {
if (arr.isEmpty()) return 0.0
val sorted = arr.copyOf().also { it.sort() }
val med = sorted[floor(sorted.size / 2.0).toInt()]
val devSorted = DoubleArray(sorted.size) { abs(sorted[it] - med) }
devSorted.sort()
return devSorted[floor(devSorted.size / 2.0).toInt()]
}
// ─────────────────────────────────────────────────────────
// Soft threshold
// ─────────────────────────────────────────────────────────
fun softThreshold(coeffs: DoubleArray, lambda: Double): DoubleArray {
val out = DoubleArray(coeffs.size)
for (i in coeffs.indices) {
val v = coeffs[i]
out[i] = when {
v > lambda -> v - lambda
v < -lambda -> v + lambda
else -> 0.0
}
}
return out
}
// ─────────────────────────────────────────────────────────
// Sun 2024 energy-ratio adaptive denoise
// ─────────────────────────────────────────────────────────
fun denoise(signal: DoubleArray, levels: Int = 3, eps: Double = 0.1): DoubleArray {
val n0 = signal.size
if (n0 < (1 shl levels)) return signal.copyOf()
val decomp = dwt(signal, levels)
val sigmas = DoubleArray(decomp.details.size) { mad(decomp.details[it]) * 1.4826 }
val energies = DoubleArray(decomp.details.size) { idx ->
val d = decomp.details[idx]
var s = 0.0
for (v in d) s += v * v
s / d.size
}
var eMax = 1e-30
for (e in energies) if (e > eMax) eMax = e
val ratios = DoubleArray(energies.size) { energies[it] / eMax }
val denoisedDetails = decomp.details.mapIndexed { j, d ->
val nj = d.size
val lambda = sigmas[j] *
sqrt(2.0 * ln(max(2.0, nj.toDouble()))) *
((1.0 - ratios[j]) + eps)
softThreshold(d, lambda)
}
val newDecomp = Decomposition(
approx = decomp.approx,
details = denoisedDetails,
lengths = decomp.lengths,
paddedLength = decomp.paddedLength
)
val reconstructed = idwt(newDecomp)
return DoubleArray(n0) { reconstructed[it] }
}
// ─────────────────────────────────────────────────────────
// Per-band reconstruction (frequency-decomposition view)
// ─────────────────────────────────────────────────────────
/**
* @param levelIdx 0..levels-1 → reconstruct only that detail level
* -1 → reconstruct only the approximation
*/
fun reconstructLevel(decomp: Decomposition, levelIdx: Int, originalLength: Int = -1): DoubleArray {
val zerosA = DoubleArray(decomp.approx.size)
val zerosD = decomp.details.map { DoubleArray(it.size) }
val stub = Decomposition(
approx = if (levelIdx == -1) decomp.approx else zerosA,
details = decomp.details.mapIndexed { j, d -> if (j == levelIdx) d else zerosD[j] },
lengths = decomp.lengths,
paddedLength = decomp.paddedLength
)
val recon = idwt(stub)
val n = if (originalLength >= 0) originalLength else recon.size
return DoubleArray(n) { recon[it] }
}
/**
* Boundary-clean per-band reconstruction with symmetric reflection padding.
* Suppresses periodic-boundary wraparound artefact at samples ~85..99.
*/
fun reconstructLevelClean(signal: DoubleArray, levels: Int, levelIdx: Int): DoubleArray {
val n = signal.size
val minPad = 32
val mLevel = 1 shl levels
val mTarget = (ceil((n + 2 * minPad).toDouble() / mLevel) * mLevel).toInt()
val padLeft = ((mTarget - n) / 2)
val padRight = mTarget - n - padLeft
val ext = DoubleArray(mTarget)
for (i in 0 until padLeft) ext[i] = signal[min(padLeft - 1 - i, n - 1)]
for (i in 0 until n) ext[padLeft + i] = signal[i]
for (i in 0 until padRight) ext[padLeft + n + i] = signal[max(n - 1 - i, 0)]
val decomp = dwt(ext, levels)
val zerosA = DoubleArray(decomp.approx.size)
val zerosD = decomp.details.map { DoubleArray(it.size) }
val stub = Decomposition(
approx = if (levelIdx == -1) decomp.approx else zerosA,
details = decomp.details.mapIndexed { j, d -> if (j == levelIdx) d else zerosD[j] },
lengths = decomp.lengths,
paddedLength = decomp.paddedLength
)
val recon = idwt(stub)
return DoubleArray(n) { recon[padLeft + it] }
}
// ─────────────────────────────────────────────────────────
// Diagnostic (no thresholding — for visualisation)
// ─────────────────────────────────────────────────────────
fun diagnose(signal: DoubleArray, levels: Int = 3): Diagnose? {
if (signal.size < (1 shl levels)) return null
val decomp = dwt(signal, levels)
val sigmas = DoubleArray(decomp.details.size) { mad(decomp.details[it]) * 1.4826 }
val energies = DoubleArray(decomp.details.size) { idx ->
val d = decomp.details[idx]
var s = 0.0
for (v in d) s += v * v
s / d.size
}
var eMax = 1e-30
for (e in energies) if (e > eMax) eMax = e
val ratios = DoubleArray(energies.size) { energies[it] / eMax }
return Diagnose(
levels = levels,
sigma = sigmas,
energy = energies,
ratio = ratios,
detailLengths = IntArray(decomp.details.size) { decomp.details[it].size }
)
}
}
@@ -0,0 +1,45 @@
/*
* Method D config — port of vesiscan_test/library/method_d/config_d.py.
*
* 모든 디폴트값을 python 과 동일하게 유지. 튜닝 근거 주석은 원본 참조.
* 각도 보정(otsu_ratio × cos(angle))은 호출부(MethodDRunner)에서 곱함.
*/
package com.medithings.vesiscan.walldetect.algo.methodd
data class MethodDParams(
val otsuRatio: Double = 0.88,
val lowMinLen: Int = 3,
val mergeGapMax: Int = 3,
val dMax: Int = 10,
val antDMax: Int = 18,
val postMaxIdx: Int = 100,
val minUrineLen: Int = 10,
val distDecay: Double = 0.1,
val promGamma: Double = 1.5,
val shoulderDistDecay: Double = 0.3,
// 02bed02 (4494098 VBTWD201): shoulder 후보 prominence 가중 1.0→1.5
// peak(promGamma=1.5) 와 정렬해 강한 shoulder 가 약한 peak 에 지지 않음.
// peak↔shoulder 재분류 토글 제거 (v2 cycle CV 12.8→9.5%).
val shoulderPromGamma: Double = 1.5,
val minShoulderProm: Double = 100.0,
val minPeakProm: Double = 50.0,
val shoulderScoreHandicap: Double = 0.15,
val gapPeakMinProm: Double = 50.0,
val minWallLumenRatio: Double = 1.16,
val minPostRawRatio: Double = 1.08,
val inwardWalkWin: Int = 3,
val inwardWalkSlopeTol: Double = 10.0,
// 02bed02 (4494098 VBTWD201): OS-CFAR window 5→9
// heavy median 강화로 span_e 안정화 → post 벽 검출 일관성 (CH0 post SD 1.45→0.36).
// urine→후벽 전이부 평탄면이 threshold 를 스쳐 span_e 가 cycle 마다 튀던 것 제거.
val oscfarWin: Int = 9,
val oscfarMaxIters: Int = 4,
val minLightSeparability: Double = 0.75,
val minSpanLen: Int = 5,
val recoveryExtendOutward: Boolean = false,
val wallRatioUsePostOnly: Boolean = true,
) {
companion object {
val DEFAULT = MethodDParams()
}
}
@@ -0,0 +1,26 @@
/*
* Method D preprocessing — port of method_d/preprocessing.py.
*
* heavy = SG (7,3) + OS-CFAR iterative median (span 검출용; ringing/speckle 흡수)
* light = SG (7,3) only (wall peak / subsample refine 용)
*
* Python config_6ch.SG_WIN=7, SG_POLY=3 — V4.1 의 (5,2) 하드코딩 커널과 분리.
* 일반화된 SgSmoothGeneric 으로 호출 (any window, polyorder 지원).
*/
package com.medithings.vesiscan.walldetect.algo.methodd
import com.medithings.vesiscan.walldetect.algo.MedianFilter
import com.medithings.vesiscan.walldetect.algo.SgSmoothGeneric
object MethodDPreprocessing {
const val SG_WIN = 7 // config_6ch.SG_WIN
const val SG_POLY = 3 // config_6ch.SG_POLY
fun preprocessHeavy(raw: DoubleArray, params: MethodDParams = MethodDParams.DEFAULT): DoubleArray {
val sg = SgSmoothGeneric.smooth(raw, SG_WIN, SG_POLY)
return MedianFilter.runningMedianRoot(sg, params.oscfarWin, params.oscfarMaxIters)
}
fun preprocessLight(raw: DoubleArray): DoubleArray =
SgSmoothGeneric.smooth(raw, SG_WIN, SG_POLY)
}
@@ -0,0 +1,107 @@
/*
* Method D span detection — port of method_d/span.py.
*
* extract_low_echo_span:
* 1) Otsu threshold × ratio = low_amp
* 2) low_mask = sg ≤ low_amp → contiguous spans (len ≥ low_min_len)
* 3) merge_close_spans (gap_peak_min_prom 가드)
* 4) 첫 candidate (sig_end 제외, len ≥ min_span_len) 채택
* 5) post 후위 검색 상한(post_max_idx) 초과 시 reject
* 6) walk_inward_to_valley 로 양쪽 valley plateau 시작점까지 shrink
* 7) walk 결과가 min_span_len 미만이면 reject
*
* extract_low_echo_span_with_fallback:
* heavy primary; heavy 가 fail 하거나 끝까지 흐르면 light 로 재시도.
* light 의 Otsu separability < min_light_separability 이면 fallback 거부.
*/
package com.medithings.vesiscan.walldetect.algo.methodd
import com.medithings.vesiscan.walldetect.algo.Otsu
import com.medithings.vesiscan.walldetect.algo.SpanUtils
object MethodDSpan {
data class SpanResult(
val lowStart: Int,
val lowEnd: Int,
val lowAmp: Double,
val inwardWalkAnt: Int,
val inwardWalkPost: Int,
)
private fun walkInwardToValley(
sg: DoubleArray, spanS: Int, spanE: Int,
win: Int, slopeTol: Double,
): IntArray {
var s = spanS
var e = spanE
val threshold = slopeTol * win
while (s + win <= e && (sg[s] - sg[s + win]) >= threshold) s++
while (e - win >= s && (sg[e] - sg[e - win]) >= threshold) e--
return intArrayOf(s, e)
}
fun extractLowEchoSpan(
sg: DoubleArray,
otsuRatio: Double,
params: MethodDParams = MethodDParams.DEFAULT,
): SpanResult? {
val otsuRes = Otsu.otsu1dWithSeparability(sg)
val lowAmp = otsuRes.threshold * otsuRatio
val mask = BooleanArray(sg.size) { sg[it] <= lowAmp }
val rawSpans = SpanUtils.contiguousTrueSpans(mask)
.filter { (it.end - it.start + 1) >= params.lowMinLen }
val spans = SpanUtils.mergeCloseSpans(
rawSpans,
maxGap = params.mergeGapMax,
sg = sg,
gapPeakThr = lowAmp + params.gapPeakMinProm,
)
if (spans.isEmpty()) return null
// 앞쪽 span 부터 순회: sig_end 제외 + 길이 ≥ min_span_len 인 첫 span.
val sigEnd = sg.size - 1
var pickStart = -1
var pickEnd = -1
for (span in spans) {
if (span.end >= sigEnd) continue
if ((span.end - span.start + 1) < params.minSpanLen) continue
pickStart = span.start
pickEnd = span.end
break
}
if (pickStart < 0) return null
if (pickEnd >= params.postMaxIdx) return null
val walked = walkInwardToValley(
sg, pickStart, pickEnd,
params.inwardWalkWin, params.inwardWalkSlopeTol,
)
val s = walked[0]
val e = walked[1]
if ((e - s + 1) < params.minSpanLen) return null
return SpanResult(
lowStart = s,
lowEnd = e,
lowAmp = lowAmp,
inwardWalkAnt = s - pickStart,
inwardWalkPost = pickEnd - e,
)
}
fun extractLowEchoSpanWithFallback(
sgHeavy: DoubleArray,
sgLight: DoubleArray,
otsuRatio: Double,
params: MethodDParams = MethodDParams.DEFAULT,
): SpanResult? {
val primary = extractLowEchoSpan(sgHeavy, otsuRatio, params)
if (primary != null && primary.lowEnd < sgHeavy.size - 1) return primary
// Light fallback — unimodal 신호 차단.
val sep = Otsu.otsu1dWithSeparability(sgLight).separability
if (sep < params.minLightSeparability) return null
return extractLowEchoSpan(sgLight, otsuRatio, params)
}
}
@@ -0,0 +1,87 @@
/*
* TGC (Time Gain Compensation) — port of denoising.py:apply_tgc_pipeline.
*
* 깊이가 깊을수록 음향 신호가 감쇠하는 현상을 보정:
* 1) fit_attenuation_lines: per-channel linear LS fit (slope, intercept) on x=[0..N-1]
* 2) adaptive_tgc_ratio(slope, slope_thresh=3.0, slope_max=15.0, ratio_min=0.1):
* |slope| < 3.0 → ratio=1.0 (보정 없음, 가파르지 않은 채널)
* else ratio = 1.0 - (1-ratio_min) * (|slope|-slope_thresh) / (slope_max-slope_thresh)
* clipped to [ratio_min, 1.0]
* 3) target_slope = slope * ratio
* 4) compensation = (target_slope - slope) * x
* 5) compensated = original + compensation
*
* slope > 0 (양수, 깊을수록 밝아짐) 인 채널은 보정 스킵 — 비정상 케이스.
*
* 입력은 (n_ch, n_samples) 단일 scan. method_d 가 heavy/light 각각에 호출.
*/
package com.medithings.vesiscan.walldetect.algo.methodd
import kotlin.math.abs
object MethodDTgc {
/** Linear LS fit y = a + b*x on x=[0..n-1] → (slope=b, intercept=a). numpy.polyfit(x,y,1). */
fun fitAttenuationLine(y: DoubleArray): Pair<Double, Double> {
val n = y.size
if (n < 2) return 0.0 to (if (n == 1) y[0] else 0.0)
// x = 0..n-1
val sumX = (n - 1).toDouble() * n / 2.0 // Σx
val sumX2 = (n - 1).toDouble() * n * (2 * n - 1) / 6.0 // Σx²
var sumY = 0.0
var sumXY = 0.0
for (i in 0 until n) {
sumY += y[i]
sumXY += i * y[i]
}
val meanX = sumX / n
val meanY = sumY / n
val varX = sumX2 - n * meanX * meanX
val covXY = sumXY - n * meanX * meanY
val slope = if (abs(varX) < 1e-12) 0.0 else covXY / varX
val intercept = meanY - slope * meanX
return slope to intercept
}
/** adaptive_tgc_ratio(slope) — Python 그대로. slope_thresh=3.0, slope_max=15.0, ratio_min=0.1. */
fun adaptiveTgcRatio(
slope: Double,
slopeThresh: Double = 3.0,
slopeMax: Double = 15.0,
ratioMin: Double = 0.1,
): Double {
val absSlope = abs(slope)
if (absSlope < slopeThresh) return 1.0
val ratio = 1.0 - (1.0 - ratioMin) * (absSlope - slopeThresh) / (slopeMax - slopeThresh)
return maxOf(ratio, ratioMin)
}
/**
* Apply TGC to a single scan (n_ch × n_samples). 채널별 slope 계산 → ratio →
* compensation = (target_slope - slope) * x 가산. slope > 0 인 채널은 skip.
*
* Python apply_tgc_pipeline(df, n_ch, center_ch=None) with center_ch=None
* defaults to all channels — 우리는 입력 list 전체에 적용.
*
* @param channels (n_ch) 길이의 (n_samples) DoubleArray
* @param ratioMin Python default 0.1. 1차 비교에선 그대로.
* @param targetRatio override (Python `target_ratio` param). null 이면 adaptive.
*/
fun applyTgcPipeline(
channels: List<DoubleArray>,
ratioMin: Double = 0.1,
targetRatio: Double? = null,
): List<DoubleArray> {
if (channels.isEmpty()) return channels
return channels.map { row ->
val (slope, _) = fitAttenuationLine(row)
if (slope >= 0.0) return@map row.copyOf()
val ratio = targetRatio ?: adaptiveTgcRatio(slope, ratioMin = ratioMin)
val targetSlope = slope * ratio
val delta = targetSlope - slope // 음수 slope → 음수 ratio 곱하면 더 작은 음수,
// delta = targetSlope - slope > 0 → 깊이 갈수록 보정+
val out = DoubleArray(row.size) { i -> row[i] + delta * i }
out
}
}
}
@@ -0,0 +1,167 @@
/*
* Method D wall selection — port of method_d/wall_select.py.
*
* span edge 근방 d_max 안에서 두 종류 후보를 수집:
* 1) local maxima (PeakDetection.findPeaks1D) → type='peak'
* 2) d2 local-min 이면서 음수 (어깨) → type='shoulder'
* ±1 sample 이내 peak 와 중복이면 제거.
*
* 스코어링:
* score = prom^gamma / (1 + dist_decay × dist)
* prom = sig[peak] - sig[adjacent_valley] (valley = peak ↔ edge 사이 최소점)
* dist = |peak - edge|
*
* shoulder 는 prom_gate(min_shoulder_prom) 더 strict,
* score 에 (1+shoulder_score_handicap) deadband 페널티 → peak 와 동률 토글 차단.
*/
package com.medithings.vesiscan.walldetect.algo.methodd
import com.medithings.vesiscan.walldetect.algo.PeakDetection
import kotlin.math.abs
import kotlin.math.pow
object MethodDWallSelect {
enum class Side { ANT, POST }
enum class CType { PEAK, SHOULDER }
data class Candidate(
val idx: Int,
val prom: Double,
val dist: Int,
val valleyIdx: Int,
val score: Double,
val type: CType,
)
data class WallResult(
val best: Candidate?,
val candidates: List<Candidate>,
)
private fun adjacentValley(sig: DoubleArray, peakIdx: Int, edgeIdx: Int): Int {
val step = if (edgeIdx > peakIdx) 1 else -1
var i = peakIdx
while (true) {
val nxt = i + step
if ((step > 0 && nxt > edgeIdx) || (step < 0 && nxt < edgeIdx)) return edgeIdx
if (sig[nxt] <= sig[i]) i = nxt
else return i
}
}
private fun findShouldersLocal(sub: DoubleArray): IntArray {
if (sub.size < 5) return IntArray(0)
val n = sub.size
val d2 = DoubleArray(n - 2) { i -> sub[i] - 2 * sub[i + 1] + sub[i + 2] }
val out = mutableListOf<Int>()
for (i in 1 until d2.size - 1) {
if (d2[i] < 0 && d2[i] < d2[i - 1] && d2[i] < d2[i + 1]) {
out += i + 1 // d2 index → sub index (+1 from central-diff offset)
}
}
return out.toIntArray()
}
fun findWallPeakLocal(
sig: DoubleArray,
spanS: Int,
spanE: Int,
side: Side,
params: MethodDParams = MethodDParams.DEFAULT,
dMaxOverride: Int? = null,
inwardWalk: Int = 0,
): WallResult {
val n = sig.size
val baseDMax = dMaxOverride ?: when (side) {
Side.ANT -> params.antDMax
Side.POST -> params.dMax
}
val effectiveDMax = baseDMax + inwardWalk
val edge: Int
val lo: Int
val hi: Int
when (side) {
Side.ANT -> {
edge = spanS
lo = (edge - effectiveDMax).coerceAtLeast(0)
hi = edge
}
Side.POST -> {
edge = spanE
lo = edge
hi = minOf(n - 1, edge + effectiveDMax, params.postMaxIdx)
}
}
if (hi <= lo) return WallResult(null, emptyList())
val sub = DoubleArray(hi - lo + 1) { sig[lo + it] }
val relPeaks = PeakDetection.findPeaks1D(sub, 0, sub.size)
val relShoulders = findShouldersLocal(sub)
val candIdxType = mutableListOf<Pair<Int, CType>>()
for (p in relPeaks) candIdxType += Pair(lo + p, CType.PEAK)
for (s in relShoulders) {
val idx = lo + s
if (candIdxType.any { it.second == CType.PEAK && abs(it.first - idx) <= 1 }) continue
candIdxType += Pair(idx, CType.SHOULDER)
}
val scored = mutableListOf<Candidate>()
for ((p, ctype) in candIdxType) {
val vIdx = adjacentValley(sig, p, edge)
val prom = sig[p] - sig[vIdx]
val gate = if (ctype == CType.PEAK) params.minPeakProm else params.minShoulderProm
if (prom <= 0 || prom < gate) continue
val dist = abs(p - edge)
val score = if (ctype == CType.SHOULDER) {
prom.pow(params.shoulderPromGamma) /
(1.0 + params.shoulderDistDecay * dist) /
(1.0 + params.shoulderScoreHandicap)
} else {
prom.pow(params.promGamma) / (1.0 + params.distDecay * dist)
}
scored += Candidate(p, prom, dist, vIdx, score, ctype)
}
if (scored.isEmpty()) return WallResult(null, emptyList())
return WallResult(scored.maxBy { it.score }, scored)
}
/**
* Parabolic sub-sample refine (Cespedes 1995) for `peak=true`.
* Returns the original idx as Double when refinement is invalid.
*/
fun refineParabolic(envelope: DoubleArray, idx: Int, peak: Boolean = true): Double {
val n = envelope.size
if (idx < 1 || idx > n - 2) return idx.toDouble()
val ym1 = envelope[idx - 1]
val y0 = envelope[idx]
val yp1 = envelope[idx + 1]
val denom = ym1 - 2 * y0 + yp1
if (denom == 0.0 || denom.isNaN() || denom.isInfinite()) return idx.toDouble()
if (peak && denom > 0) return idx.toDouble()
if (!peak && denom < 0) return idx.toDouble()
val delta = 0.5 * (ym1 - yp1) / denom
if (delta.isNaN() || delta.isInfinite() || abs(delta) > 1.0) return idx.toDouble()
return idx + delta
}
/**
* Shoulder (d2 local-min) sub-sample refine via parabolic fit on d2 itself.
* Needs envelope[idx-2 .. idx+2].
*/
fun refineShoulder(envelope: DoubleArray, idx: Int): Double {
val n = envelope.size
if (idx < 2 || idx > n - 3) return idx.toDouble()
val e = envelope
val d2m1 = e[idx - 2] - 2 * e[idx - 1] + e[idx]
val d20 = e[idx - 1] - 2 * e[idx] + e[idx + 1]
val d2p1 = e[idx] - 2 * e[idx + 1] + e[idx + 2]
val denom = d2m1 - 2 * d20 + d2p1
if (denom <= 0.0 || denom.isNaN() || denom.isInfinite()) return idx.toDouble()
val delta = 0.5 * (d2m1 - d2p1) / denom
if (delta.isNaN() || delta.isInfinite() || abs(delta) > 1.0) return idx.toDouble()
return idx + delta
}
}
@@ -0,0 +1,94 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/core/config.js (1:1).
* 상수 변경 금지 — 정의가 바뀌면 V2/V4.1 골든 테스트가 깨집니다.
*/
package com.medithings.vesiscan.walldetect.core
object WdConfig {
const val SG_WIN = 5
const val SG_POLY = 2
const val OSCFAR_WIN = 15
const val OSCFAR_RANK = 0.4
const val OSCFAR_SCALE = 1.05
const val CHORD_IDEAL_MM = 50.0
const val CHORD_SIGMA_MM = 30.0
val CHORD_MM_RANGE = doubleArrayOf(5.0, 130.0)
const val PROMINENCE_DB_MIN = 1.5
const val W_LUMEN_DARKNESS = 25.0
const val W_CHORD_PRIOR = 8.0
const val W_SYMMETRY = 3.0
const val CONTRAST_ALPHA = 2.0
const val CONTRAST_TAU = 6
const val CONTRAST_WIN = 10
// 2026-04-29 #2: aligned to study/amode_simulator.html V4 운용 확정 (DPS card).
// Theory : DPS = c_ref(1530 m/s) · Δt(2.524 μs) / 2 = 1.931 mm/sp
// Phantom : D(100.406 mm) / chord(52 samples) = 1.9309 mm/sp
// K : c_eff / c_ref = 1.000
// 결과적으로 amode_simulator 와 동일 좌표계 → 530 mL phantom 시 R≈50.20mm 회복.
// Was `const val` — promoted to mutable so settings can override it at
// runtime (UI "DPS" field syncs PiezoHW.distancePerSample → DPS_DEFAULT
// so V41Detector's Geometry / AnatomicalGate use the same value the
// 6-channel BV estimator uses).
@Volatile @JvmField var DPS_DEFAULT: Double = 1.9309
const val DELAY_MM_DEFAULT = 6.85
const val C_EFF_DEFAULT = 1530.0
// ── Low-echo detection tuning (method_b SSOT, config_6ch.py L116-127) ──
const val LOW_ECHO_AMP = 1250.0
const val LOW_MIN_LEN = 3
const val MERGE_GAP_MAX = 3
const val PEAK_SEARCH_WIN = 20
const val POST_MAX_IDX = 80
const val MIN_PEAK_MARGIN = 30.0
const val MIN_URINE_LEN = 3
const val EDGE_DIST_DECAY = 0.12
const val MAX_PEAK_CANDIDATES = 3
const val VALLEY_STOP_RISE = 50.0
const val PEAK_LO = 4
const val PEAK_HI_MARGIN = 2
const val PHANTOM_530_ANT_DEPTH_MM = 32.0
const val PHANTOM_530_R_MM = 50.20
const val PHANTOM_530_BV_ML = 530.0
const val USE_LR_TILT = false
}
object WdProbe {
// 2026-04-29 #2: aligned to study/amode_simulator.html v2 (30°) probe.
// Mechanical SI : [0, -10, -20, -30, -10, -10] (CAD)
// Snell SI : [0.0, -6.89, -13.66, -20.20, -6.87, -6.87] ← acoustic ray
// Mechanical LR : [0, 0, 0, 0, -5, 5]
// Snell LR : [0, 0, 0, 0, -3.42, 3.42]
// sensor_z (mm) : [19.3, 13.0, 6.7, 0.0, 9.85, 9.85] CH0=top, CH3=bottom
// sensor_x (mm) : [0, 0, 0, 0, -10, 10]
val DEGREE = doubleArrayOf(0.0, -6.89, -13.66, -20.20, -6.87, -6.87)
val DEGREE_LR = doubleArrayOf(0.0, 0.0, 0.0, 0.0, -3.42, 3.42)
val MECH_DEGREE = doubleArrayOf(0.0, -10.0, -20.0, -30.0, -10.0, -10.0)
val MECH_DEGREE_LR = doubleArrayOf(0.0, 0.0, 0.0, 0.0, -5.0, 5.0)
val SENSOR_Z = doubleArrayOf(19.3, 13.0, 6.7, 0.0, 9.85, 9.85)
val SENSOR_X = doubleArrayOf(0.0, 0.0, 0.0, 0.0, -10.0, 10.0)
/** v1 — 10° device (alternate). */
val V1_DEGREE = doubleArrayOf(6.89, 0.0, -6.89, -13.66, -6.87, -6.87)
val V1_DEGREE_LR = doubleArrayOf(0.0, 0.0, 0.0, 0.0, -3.42, 3.42)
val V1_SENSOR_Z = doubleArrayOf(20.1, 12.8, 7.0, 0.0, 9.9, 9.9)
val PIEZO_W = intArrayOf(12, 12, 12, 12, 6, 6)
val PIEZO_H = intArrayOf(6, 6, 6, 6, 12, 12)
const val PIEZO_T = 1.0
val Z_RECESS = doubleArrayOf(2.37, 1.188, 1.188, 1.188, 0.0, 0.0)
const val HOUSING_W = 30
const val HOUSING_H = 40
const val HOUSING_DEPTH = 8
const val LAYOUT_NOTE = "CH0 (top) → CH3 (navel) vertical column · CH4 = left · CH5 = right"
}
@@ -0,0 +1,98 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/core/linalg.js (1:1).
* Small dense linear-algebra helpers (Gauss-Jordan solve with partial pivoting).
* Used by sphere_kasa / sphere_lm.
*/
package com.medithings.vesiscan.walldetect.core
import kotlin.math.abs
import kotlin.math.sqrt
object WdLinalg {
/**
* Solve M·x = b where M is n×n. Returns x or null on singular matrix.
* Mutates internal copies; original M and b are untouched.
*/
fun solve(m: Array<DoubleArray>, b: DoubleArray): DoubleArray? {
val n = b.size
// Deep copy
val a = Array(m.size) { m[it].copyOf() }
val x = b.copyOf()
for (i in 0 until n) {
// Partial pivot
var p = i
for (r in (i + 1) until n) {
if (abs(a[r][i]) > abs(a[p][i])) p = r
}
if (abs(a[p][i]) < 1e-12) return null
if (p != i) {
val tmpRow = a[i]; a[i] = a[p]; a[p] = tmpRow
val tmp = x[i]; x[i] = x[p]; x[p] = tmp
}
val piv = a[i][i]
for (j in i until n) a[i][j] /= piv
x[i] /= piv
for (r in 0 until n) {
if (r == i) continue
val f = a[r][i]
if (f == 0.0) continue
for (j in i until n) a[r][j] -= f * a[i][j]
x[r] -= f * x[i]
}
}
return x
}
/** Identity n×n. */
fun eye(n: Int): Array<DoubleArray> {
val ii = Array(n) { DoubleArray(n) }
for (i in 0 until n) ii[i][i] = 1.0
return ii
}
/**
* out = A^T · A.
* A is m×n (rows of length n).
*/
fun ata(a: Array<DoubleArray>): Array<DoubleArray> {
val m = a.size
val n = a[0].size
val out = Array(n) { DoubleArray(n) }
for (i in 0 until n) {
for (j in 0 until n) {
var s = 0.0
for (k in 0 until m) s += a[k][i] * a[k][j]
out[i][j] = s
}
}
return out
}
/** A^T · b (m×n × m → n). */
fun atb(a: Array<DoubleArray>, b: DoubleArray): DoubleArray {
val m = a.size
val n = a[0].size
val out = DoubleArray(n)
for (i in 0 until n) {
var s = 0.0
for (k in 0 until m) s += a[k][i] * b[k]
out[i] = s
}
return out
}
fun vecNorm(v: DoubleArray): Double {
var s = 0.0
for (x in v) s += x * x
return sqrt(s)
}
}
@@ -0,0 +1,63 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Mirror of study/wall_detect_verify/js/core/numeric.js (1:1).
* Pure numeric helpers (median / quantile / mean / std / minmax).
*
* 정밀도: 내부 계산은 Double로 수행 (JS는 모두 double). 호출자가 Float 으로 다운캐스트.
*/
package com.medithings.vesiscan.walldetect.core
import kotlin.math.ceil
import kotlin.math.floor
import kotlin.math.sqrt
object WdNumeric {
/**
* Linear-interpolated quantile, matching JS `quantile()`:
* pos = q * (n - 1); lo = floor(pos); hi = ceil(pos);
* a[lo] + (a[hi] - a[lo]) * (pos - lo)
*/
fun quantile(arr: DoubleArray, q: Double): Double {
if (arr.isEmpty()) return Double.NaN
val a = arr.copyOf().also { it.sort() }
val pos = q * (a.size - 1)
val lo = floor(pos).toInt()
val hi = ceil(pos).toInt()
if (lo == hi) return a[lo]
return a[lo] + (a[hi] - a[lo]) * (pos - lo)
}
fun median(arr: DoubleArray): Double = quantile(arr, 0.5)
fun mean(arr: DoubleArray): Double {
if (arr.isEmpty()) return Double.NaN
var s = 0.0
for (v in arr) s += v
return s / arr.size
}
fun std(arr: DoubleArray, m: Double? = null): Double {
if (arr.isEmpty()) return Double.NaN
val mu = m ?: mean(arr)
var s = 0.0
for (v in arr) {
val d = v - mu
s += d * d
}
return sqrt(s / arr.size)
}
fun minmax(arr: DoubleArray): DoubleArray {
var lo = Double.POSITIVE_INFINITY
var hi = Double.NEGATIVE_INFINITY
for (v in arr) {
if (v < lo) lo = v
if (v > hi) hi = v
}
return doubleArrayOf(lo, hi)
}
}
@@ -0,0 +1,38 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Multi-channel BV estimation result (4-layer dispatch — see
* docs/BV-CALCULATION-DESIGN.md). Replaces the sphere-only `v41Sphere.bvMl`
* for the primary "BV" metric while sphere fit stays for cross-check.
*/
package com.medithings.vesiscan.walldetect.dto
data class BvDispatchResult(
val bvMl: Float?, // null = no estimate
val rMm: Float?, // equivalent sphere radius (sphere or chord-derived)
val method: String, // FrustumLR | FrustumNoLR | ChordMedian | Verathon | None
val confidence: Float, // 0..1
val nCenter: Int, // gated center channels (CH0..CH3) AFTER consensus filter
val nLateral: Int, // gated lateral channels (CH4, CH5) AFTER consensus filter
val lrRatio: Float, // applied LR/AP ratio (1.0 if no lateral)
val warnings: List<String> = emptyList(),
val sphereCrossCheckBvMl: Float? = null,
// ── ChordConsensus integration (v4.1.1) ──────────────────────────────
/** Trusted channel indices after score ≥ 0.40 + Tukey/Fischler-Bolles
* consensus filter. Empty when no detections passed. */
val trustedChannels: List<Int> = emptyList(),
/** Channels rejected by the consensus filter, with reason string. */
val rejectedChannels: List<RejectedChannel> = emptyList(),
/** Median chord across the trusted set (mm). */
val consensusMedianChordMm: Float? = null,
/** MAD of trusted chords (mm). null when N < 3. */
val consensusMadMm: Float? = null,
/** Highest-score trusted channel (RANSAC leader). */
val leaderCh: Int? = null,
)
data class RejectedChannel(val ch: Int, val chordMm: Float, val reason: String)
@@ -0,0 +1,27 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Per-channel detector result. Schema is identical for V2 and V4.1.
* V2 결과는 v41Diag 가 항상 null.
*/
package com.medithings.vesiscan.walldetect.dto
data class ChannelResult(
val ch: Int,
val sg: List<Float>,
val threshold: List<Float>,
val antIdx: Int?,
val postIdx: Int?,
val antRefined: Float?,
val postRefined: Float?,
val antMm: Float?,
val postMm: Float?,
val lumenStart: Int?,
val lumenEnd: Int?,
val chordMm: Float?,
val clipping: Boolean? = null,
val v41Diag: V41Diagnostics? = null
)
@@ -0,0 +1,20 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Top-level detector output. Identical schema for V2 and V4.1.
* Only `algorithm` / `algorithmVersion` differ; all other field names match.
*/
package com.medithings.vesiscan.walldetect.dto
data class DetectionResult(
val requestId: String,
val timestampMs: Long,
val algorithm: String,
val algorithmVersion: String,
val processingMs: Double,
val perChannel: List<ChannelResult>,
val summary: DetectionSummary
)
@@ -0,0 +1,24 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Sweep-level summary. Schema identical for V2 and V4.1.
* V4.1-only fields are null in V2 results.
*/
package com.medithings.vesiscan.walldetect.dto
data class DetectionSummary(
val matchCount: Int,
val chordMmMean: Float?,
val chordMmStd: Float?,
val tier: String? = null,
val scoreMean: Float? = null,
val scoreMin: Float? = null,
val scoreMax: Float? = null,
val gatedCount: Int? = null,
val v41Sphere: V41SphereFit? = null,
/** Multi-channel BV dispatch (PR-13). Replaces sphere-only BV as primary value. */
val bvDispatch: BvDispatchResult? = null,
)
@@ -0,0 +1,29 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* Common input DTO for both V2 and V4.1 detectors.
* 두 detector 가 동일 입력을 받는다는 drift-zero 보장의 단일 source.
*/
package com.medithings.vesiscan.walldetect.dto
data class SweepInput(
val requestId: String,
val timestampMs: Long,
val deviceId: String,
val sweepSeq: Int,
val probe: ProbeProfileDto,
val adc: List<List<Int>>,
val gainDb: Int = 0
)
data class ProbeProfileDto(
val name: String,
val rMm: Double,
val anteriorAnchorMm: Double,
val fsHz: Long,
val cMmPerUs: Double,
val anglesDeg: List<Double>
)
@@ -0,0 +1,41 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* V4.1-only per-channel diagnostics carrier.
* V2 결과에서는 ChannelResult.v41Diag = null. 스키마 자체는 양쪽 동일.
*/
package com.medithings.vesiscan.walldetect.dto
data class V41Diagnostics(
val waveletDenoised: List<Float>,
val waveletEnergyRatio: WaveletEnergyRatio,
val waveletCoefs: WaveletCoefs,
val peaksAll: List<Int>,
val spans: List<IntRange2>,
val score: Float,
val scoreSub: ScoreSubscores,
val gated: Boolean,
val tier: String,
val sContrast: Float,
val sContrastTier: String,
val singleChannelBvMl: Float?
)
data class WaveletEnergyRatio(val l1: Float, val l2: Float, val l3: Float)
data class WaveletCoefs(
val l1: List<Float>,
val l2: List<Float>,
val l3: List<Float>,
val a3: List<Float>
)
data class ScoreSubscores(
val uFarPost: Float,
val uLumDark: Float,
val uAntGrad: Float,
val uPostGrad: Float
)
@@ -0,0 +1,22 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*
* V4.1-only sweep-level sphere fit + BV result.
* Q5: gated < 4 인 sweep 에서는 null (JS sphere_fit_2step.js 와 동일).
*/
package com.medithings.vesiscan.walldetect.dto
data class V41SphereFit(
val mode: String,
val center: Vec3,
val radiusMm: Float,
val bvMl: Float,
val residualStdMm: Float,
val wallPoints: List<Vec3>,
val deltaRMm: Float? = null,
val deltaBvMl: Float? = null,
val nPoints: Int
)
@@ -0,0 +1,11 @@
/*
* Copyright (c) 2026 Medithings Co., Ltd.
* Author: Charles KWON OhJun <charleskwon@medithings.co.kr>
* Project: CharlesKWONsLaw — wall-detect live compare
*/
package com.medithings.vesiscan.walldetect.dto
data class Vec3(val x: Float, val y: Float, val z: Float)
data class IntRange2(val start: Int, val end: Int)