feat: Method C (V4.1) 대규모 업데이트 — 2-pass pipeline + ChordConsensus

새 알고리즘 모듈 (7개):
- ImpulseReject: 2-pass Hampel filter (lumen 내 impulse 제거)
- MedianFilter: running-median pre-filter (window=7, speckle 제거)
- MorphClose: morphological closing on lumen mask (disabled)
- StaLta: STA/LTA impulse detector (Allen 1978, wall edge vs reverberation)
- ChordConsensus: Tukey MAD outlier + Fischler-Bolles pair test
- DetectionSanity: detection sanity checks
- SweepStabilizer: sweep temporal stability

변경된 모듈 (9개):
- V41Detector: 2-pass pipeline (ImpulseReject.detectWithLumenClean),
  BModeScore V41_WEIGHTS (wamp+stalta), ChordConsensus filtering
- BvEstimation: ChordConsensus filter → trusted channel only
- DetectLumenFirst: running-median pre-filter, V41_PARAMS,
  mergeGapMax=5, Kremkau half-amplitude ant fallback
- AnatomicalGate: antDepthMin 22→10mm, strict threshold 제거
- BModeScore: V41_WEIGHTS (far=0.15, dark=0.10, wamp=0.45, stalta=0.20)
- WallSelect: MAX_PEAK_CANDIDATES_POST=64
- SpanUtils, Otsu, BvDispatchResult 업데이트

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
This commit is contained in:
2026-05-04 17:01:55 +09:00
parent d6e981e1d2
commit 82634c2107
16 changed files with 1471 additions and 73 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.example.medilightv2android.walldetect
import com.example.medilightv2android.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,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.example.medilightv2android.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() }
}
@@ -28,7 +28,7 @@
* ▶ OUTPUT · `dto.DetectionResult` * ▶ OUTPUT · `dto.DetectionResult`
* ─────── * ───────
* algorithm = "v4_1" * algorithm = "v4_1"
* algorithmVersion = "v4.1.0" * algorithmVersion = "v4.1.1"
* processingMs : Double * processingMs : Double
* perChannel[6] : ChannelResult * perChannel[6] : ChannelResult
* ├─ sg[100] sg-smoothed envelope * ├─ sg[100] sg-smoothed envelope
@@ -104,10 +104,12 @@ import com.example.medilightv2android.walldetect.algo.AnatomicalGate
import com.example.medilightv2android.walldetect.algo.BModeScore import com.example.medilightv2android.walldetect.algo.BModeScore
import com.example.medilightv2android.walldetect.algo.BvEstimation import com.example.medilightv2android.walldetect.algo.BvEstimation
import com.example.medilightv2android.walldetect.algo.BvFromSphere import com.example.medilightv2android.walldetect.algo.BvFromSphere
import com.example.medilightv2android.walldetect.algo.ChordConsensus
import com.example.medilightv2android.walldetect.algo.Clipping import com.example.medilightv2android.walldetect.algo.Clipping
import com.example.medilightv2android.walldetect.algo.ContrastAux import com.example.medilightv2android.walldetect.algo.ContrastAux
import com.example.medilightv2android.walldetect.algo.DetectLumenFirst import com.example.medilightv2android.walldetect.algo.DetectLumenFirst
import com.example.medilightv2android.walldetect.algo.Geometry import com.example.medilightv2android.walldetect.algo.Geometry
import com.example.medilightv2android.walldetect.algo.ImpulseReject
import com.example.medilightv2android.walldetect.algo.PeakDetection import com.example.medilightv2android.walldetect.algo.PeakDetection
import com.example.medilightv2android.walldetect.algo.SphereFit2Step import com.example.medilightv2android.walldetect.algo.SphereFit2Step
import com.example.medilightv2android.walldetect.algo.SubsampleRefine import com.example.medilightv2android.walldetect.algo.SubsampleRefine
@@ -132,7 +134,19 @@ class V41Detector(
) : WallDetector { ) : WallDetector {
override val algorithmId: String = DetectorIds.V4_1 override val algorithmId: String = DetectorIds.V4_1
override val algorithmVersion: String = "v4.1.0" // 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 { override fun detect(input: SweepInput): DetectionResult {
val t0 = System.nanoTime() val t0 = System.nanoTime()
@@ -161,10 +175,16 @@ class V41Detector(
// Use the gated ant/post pairs (post anatomical gate) as input. Sphere // Use the gated ant/post pairs (post anatomical gate) as input. Sphere
// fit BV is passed as cross-check. // fit BV is passed as cross-check.
val bvDispatch = run { 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 dets = perCh.map { cr ->
val isGated = cr.v41Diag?.gated == true
BvEstimation.Detection( BvEstimation.Detection(
ant = if (cr.v41Diag?.gated == true) cr.antIdx else null, ant = if (isGated) cr.antIdx else null,
post = if (cr.v41Diag?.gated == true) cr.postIdx else null, post = if (isGated) cr.postIdx else null,
score = if (isGated) (cr.v41Diag?.score?.toDouble() ?: 0.0) else 0.0,
) )
} }
BvEstimation.estimate( BvEstimation.estimate(
@@ -217,8 +237,19 @@ class V41Detector(
val wlDiag = WaveletDenoise.diagnose(raw, levels = 3) val wlDiag = WaveletDenoise.diagnose(raw, levels = 3)
val wlDecomp = WaveletDenoise.dwt(raw, levels = 3) val wlDecomp = WaveletDenoise.dwt(raw, levels = 3)
// 1.a + 2.a + 2.d + 3.a + 3.b (DetectLumenFirst computes sg + per-sample CFAR + spans + walls) // 1.a + 2.a + 2.d + 3.a + 3.b — V4.1 two-pass detection.
val det = DetectLumenFirst.detect(rawInt, DetectLumenFirst.Mode.Adaptive) // • 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 sg = det.sg
val cfarThr = det.cfarThr ?: DoubleArray(sg.size) { det.adaptiveT } // safety val cfarThr = det.cfarThr ?: DoubleArray(sg.size) { det.adaptiveT } // safety
@@ -228,9 +259,14 @@ class V41Detector(
// 2.c peaks (sg, all candidates) // 2.c peaks (sg, all candidates)
val peaksAll = PeakDetection.findPeaks1D(sg).toList() val peaksAll = PeakDetection.findPeaks1D(sg).toList()
// 3.c subsample refine (parabolic, peak kind) // 3.c subsample refine (parabolic, peak kind) — operates on RAW envelope.
val antRefined = det.ant?.let { SubsampleRefine.refineParabolic(sg, it, SubsampleRefine.Kind.PEAK) } // SG smoothing slightly biases the parabola vertex, so refine uses the
val postRefined = det.post?.let { SubsampleRefine.refineParabolic(sg, it, SubsampleRefine.Kind.PEAK) } // 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() } val antMmRaw: Float? = antRefined?.let { Geometry.sampleToMm(it).toFloat() }
?: det.ant?.let { Geometry.sampleToMm(it.toDouble()).toFloat() } ?: det.ant?.let { Geometry.sampleToMm(it.toDouble()).toFloat() }
@@ -242,8 +278,17 @@ class V41Detector(
val sContrast: Float = (cR?.contrast ?: 0.0).toFloat() val sContrast: Float = (cR?.contrast ?: 0.0).toFloat()
val sContrastTier: String = ContrastAux.tierBand(cR?.contrast) val sContrastTier: String = ContrastAux.tierBand(cR?.contrast)
// 4.b B-mode composite score (RAW envelope) // 4.b B-mode composite score (RAW envelope) — V4.1 weighting.
val sR = BModeScore.score(raw, det.ant, det.post) // 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 score: Float = (sR?.total ?: 0.0).toFloat()
val scoreSub: ScoreSubscores = if (sR != null) ScoreSubscores( val scoreSub: ScoreSubscores = if (sR != null) ScoreSubscores(
uFarPost = sR.sub.far.toFloat(), uFarPost = sR.sub.far.toFloat(),
@@ -35,12 +35,36 @@ object AnatomicalGate {
val chordSlack: Double val chordSlack: Double
) )
/** Default for the 530 mL BP2 phantom (asymmetric band, admits corner captures). */ /**
val PHANTOM_530 = Preset( * 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, antDepthMin = 22.0, antDepthMax = 60.0,
rMax = 50.20, chordMin = 5.0, chordSlack = 10.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. */ /** Free-bladder clinical preset. */
val CLINICAL = Preset( val CLINICAL = Preset(
antDepthMin = 20.0, antDepthMax = 70.0, antDepthMin = 20.0, antDepthMax = 70.0,
@@ -41,13 +41,52 @@ object BModeScore {
// Defaults — DO NOT CHANGE without bumping algorithmVersion (golden tests will fail). // Defaults — DO NOT CHANGE without bumping algorithmVersion (golden tests will fail).
const val DEFAULT_TAU = 6 const val DEFAULT_TAU = 6
const val DEFAULT_WIN = 10 const val DEFAULT_WIN = 10
val DEFAULT_X50 = X50(far = 100.0, dark = 150.0, grad = 250.0) 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) val DEFAULT_WEIGHTS = Weights(far = 0.45, dark = 0.25, ant = 0.15, post = 0.15, wamp = 0.0, stalta = 0.0)
data class X50(val far: Double, val dark: Double, val grad: Double) /**
data class Weights(val far: Double, val dark: Double, val ant: Double, val post: Double) * 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 Subscores(val far: Double, val dark: Double, val antGrad: Double, val postGrad: Double) 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( data class RawValues(
val sFar: Double, val sFar: Double,
@@ -56,10 +95,16 @@ object BModeScore {
val gPost: Double, val gPost: Double,
val lumenMean: Double, val lumenMean: Double,
val farMean: Double?, val farMean: Double?,
val outsideMean: 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) 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 Tier(val tier: String, val desc: String)
@@ -131,6 +176,21 @@ object BModeScore {
} }
val outsideMean = if (oN > 0) oSum / oN else 0.0 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) // Raw subscore values (ADC)
val sFar = if (farMean != null) maxOf(0.0, farMean - lumenMean) else 0.0 val sFar = if (farMean != null) maxOf(0.0, farMean - lumenMean) else 0.0
val sDark = maxOf(0.0, outsideMean - lumenMean) val sDark = maxOf(0.0, outsideMean - lumenMean)
@@ -138,26 +198,33 @@ object BModeScore {
abs(envelope[ant + 1] - envelope[ant - 1]) / 2.0 else 0.0 abs(envelope[ant + 1] - envelope[ant - 1]) / 2.0 else 0.0
val gPost = if (post >= 1 && post <= n - 2) val gPost = if (post >= 1 && post <= n - 2)
abs(envelope[post + 1] - envelope[post - 1]) / 2.0 else 0.0 abs(envelope[post + 1] - envelope[post - 1]) / 2.0 else 0.0
val sWamp = maxOf(0.0, pMax - lumenMean)
// Mapped subscores ∈ [0,1] // Mapped subscores ∈ [0,1]
val uFar = softSat(sFar, x50.far) val uFar = softSat(sFar, x50.far)
val uDark = softSat(sDark, x50.dark) val uDark = softSat(sDark, x50.dark)
val uAnt = softSat(gAnt, x50.grad) val uAnt = softSat(gAnt, x50.grad)
val uPost = softSat(gPost, 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 + val total = weights.far * uFar +
weights.dark * uDark + weights.dark * uDark +
weights.ant * uAnt + weights.ant * uAnt +
weights.post * uPost weights.post * uPost +
weights.wamp * uWamp +
weights.stalta * uStl
val cls = classify(total) val cls = classify(total)
return Result( return Result(
total = total, total = total,
tier = cls.tier, tier = cls.tier,
desc = cls.desc, desc = cls.desc,
sub = Subscores(uFar, uDark, uAnt, uPost), sub = Subscores(uFar, uDark, uAnt, uPost, uWamp, uStl),
raw = RawValues(sFar, sDark, gAnt, gPost, lumenMean, farMean, outsideMean), raw = RawValues(sFar, sDark, gAnt, gPost, lumenMean, farMean, outsideMean,
windows = Windows(lLo = ant + 1, lHi = post, fLo = fLo, fHi = fHi) sWamp = sWamp, rStl = rStl, pMax = pMax),
windows = Windows(lLo = ant + 1, lHi = post, fLo = fLo, fHi = fHi,
pLo = pLo, pHi = pHi)
) )
} }
} }
@@ -37,12 +37,24 @@ object BvEstimation {
private const val LR_NO_DETECTION = 1.0f private const val LR_NO_DETECTION = 1.0f
private const val AREA_K = PI / 4.0 // (π/4) D² private const val AREA_K = PI / 4.0 // (π/4) D²
/** Per-channel detection input — only the integer ant/post matter here. */ /** Per-channel detection input — ant/post indices + the v4.1 B-mode score. */
data class Detection(val ant: Int?, val post: Int?) data class Detection(val ant: Int?, val post: Int?, val score: Double = 0.0)
/** /**
* Main entry point. `walls` size must be 6. * Main entry point. `walls` size must be 6. v4.1.1 dispatch:
* Returns a `BvDispatchResult`; never null (always has a method, even "None"). *
* 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( fun estimate(
walls: List<Detection>, walls: List<Detection>,
@@ -51,37 +63,105 @@ object BvEstimation {
delayMm: Double = WdConfig.DELAY_MM_DEFAULT, delayMm: Double = WdConfig.DELAY_MM_DEFAULT,
siDeg: DoubleArray = WdProbe.DEGREE, siDeg: DoubleArray = WdProbe.DEGREE,
sensorZ: DoubleArray = WdProbe.SENSOR_Z, sensorZ: DoubleArray = WdProbe.SENSOR_Z,
trustScore: Double = ChordConsensus.DEFAULT_TRUST_SCORE,
chordTol: Double = ChordConsensus.DEFAULT_CHORD_TOL,
kMad: Double = ChordConsensus.DEFAULT_K_MAD,
): BvDispatchResult { ): BvDispatchResult {
require(walls.size == 6) { "walls must have 6 entries (one per channel)" } require(walls.size == 6) { "walls must have 6 entries (one per channel)" }
val centerIdx = (0..3).filter { walls[it].ant != null && walls[it].post != null } // ── 1) ChordConsensus filter (score-trust + Tukey/Fischler-Bolles) ──
val lateralIdx = (4..5).filter { walls[it].ant != null && walls[it].post != null } 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 nC = centerIdx.size
val nL = lateralIdx.size val nL = lateralIdx.size
val nWallPts = (nC + nL) * 2
// ── Compute LR ratio (uses lateral channels if available) ── // ── LR ratio uses ONLY trusted channels too ──
val lrRatio = computeLrRatio(walls, centerIdx, lateralIdx, dps, delayMm, siDeg) val lrRatio = computeLrRatio(walls, centerIdx, lateralIdx, dps, delayMm, siDeg)
// ── Dispatch (2026-04-29 #3 — py2 _bv_core SSOT, lateral=auxiliary) ─ // ── 3) Dispatch ──
// CH4/5 (lateral) are AUXILIARY indicators only:
// • They feed `compute_lr_ratio` (S = π/4 · D² · lr_ratio)
// • They do NOT participate in the SI-axis frustum integration.
// BV is therefore always computed from the CENTER channels (CH0-3):
// • nC ≥ 2 → frustum + caps (length ∝ nC)
// • nC = 1 → Verathon chord-as-D × √lr_prior (single SI chord)
// • nC = 0 → None (no SI information for integration)
val primary = when { val primary = when {
nC == 0 -> none(nC, nL, lrRatio, nC >= 2 -> frustum(walls, centerIdx, dps, delayMm,
if (nL == 0) "no detection" else "lateral only — no SI integration") 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], nC == 1 -> verathon(walls[centerIdx[0]], centerIdx[0],
dps, siDeg, lrRatio, nC, nL, sphereCrossCheckBvMl) dps, siDeg, lrRatio, nC, nL, sphereCrossCheckBvMl)
else -> frustum(walls, centerIdx, dps, delayMm, (nC + nL) >= 2 -> chordMedian(consensus, nC, nL, lrRatio,
siDeg, sensorZ, nC, nL, lrRatio, sphereCrossCheckBvMl, "lateral-only consensus")
sphereCrossCheckBvMl) else -> none(nC, nL, lrRatio,
if ((nC + nL) == 0) "no trusted detection (consensus filter)"
else "single lateral — no SI / chord-median basis")
} }
return applyAnatomicalBounds(primary) // ── 4) Attach consensus diagnostics to the result ──
val withConsensus = primary.copy(
trustedChannels = consensus.trusted.sorted(),
rejectedChannels = consensus.rejected.map {
com.example.medilightv2android.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,
)
} }
// ────────────────────────────────────────────────────────── // ──────────────────────────────────────────────────────────
@@ -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.example.medilightv2android.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)
}
}
@@ -17,6 +17,11 @@ import com.example.medilightv2android.walldetect.core.WdNumeric
object DetectLumenFirst { 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 { object Params {
// Aligned to py2/for_app_share/config_6ch.py — LOW_ECHO_AMP bumped 1150 → 1250. // 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_ECHO_AMP = 1250.0
@@ -24,13 +29,109 @@ object DetectLumenFirst {
const val MERGE_GAP_MAX = 3 const val MERGE_GAP_MAX = 3
const val PEAK_SEARCH_WIN = 20 const val PEAK_SEARCH_WIN = 20
const val POST_MAX_IDX = 80 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_PEAK_MARGIN = 30.0
const val MIN_URINE_LEN = 3 const val MIN_URINE_LEN = 3
const val GAP_PEAK_MARGIN = 5.0 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 EDGE_DIST_DECAY = 0.12
const val VALLEY_STOP_RISE = 50.0 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 { sealed interface Mode {
/** Fixed amplitude threshold (legacy V2). */ /** Fixed amplitude threshold (legacy V2). */
data class Fixed(val thr: Double = Params.LOW_ECHO_AMP) : Mode data class Fixed(val thr: Double = Params.LOW_ECHO_AMP) : Mode
@@ -69,26 +170,46 @@ object DetectLumenFirst {
fun detect(rawAdc: IntArray, mode: Mode = Mode.Fixed()): Result { fun detect(rawAdc: IntArray, mode: Mode = Mode.Fixed()): Result {
val raw = DoubleArray(rawAdc.size) { rawAdc[it].toDouble() } val raw = DoubleArray(rawAdc.size) { rawAdc[it].toDouble() }
return detect(raw, mode) return detect(raw, mode, DEFAULT_PARAMS)
} }
fun detect(rawAdc: DoubleArray, mode: Mode = Mode.Fixed()): Result { 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 raw = rawAdc.copyOf()
val sg = Denoising.sgSmooth(raw) val sg = Denoising.sgSmooth(raw)
val n = sg.size 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 cfarThr: DoubleArray?
val T: Double val T: Double
when (mode) { when (mode) {
is Mode.Adaptive -> { is Mode.Adaptive -> {
// V4.1 OS-CFAR (Rohling 1983). // V4.1 OS-CFAR (Rohling 1983).
cfarThr = ThresholdOsCfar.perSample(sg) cfarThr = ThresholdOsCfar.perSample(sgForMask)
T = WdNumeric.median(cfarThr) T = WdNumeric.median(cfarThr)
} }
is Mode.Otsu -> { is Mode.Otsu -> {
// py2 method_b — Otsu over sg[0..POST_MAX_IDX] (skip far-tail). // py2 method_b — Otsu over sg[0..POST_MAX_IDX] (skip far-tail).
cfarThr = null cfarThr = null
val cap = minOf(Params.POST_MAX_IDX + 1, sg.size) val cap = minOf(params.postMaxIdx + 1, sg.size)
val view = DoubleArray(cap) { sg[it] } val view = DoubleArray(cap) { sg[it] }
T = Otsu.otsu1d(view) T = Otsu.otsu1d(view)
} }
@@ -102,12 +223,47 @@ object DetectLumenFirst {
} }
} }
val lowMask = BooleanArray(n) { sg[it] <= T } val lowMask = BooleanArray(n) { sgForMask[it] <= T }
val rawSpans = SpanUtils.contiguousTrueSpans(lowMask) val rawSpans = SpanUtils.contiguousTrueSpans(lowMask)
.filter { (it.end - it.start + 1) >= Params.LOW_MIN_LEN } .filter { (it.end - it.start + 1) >= params.lowMinLen }
val gapPeakThr = T + Params.GAP_PEAK_MARGIN // Stage 1 — width-bounded amplitude-aware merge.
val spans = SpanUtils.mergeCloseSpans(rawSpans, Params.MERGE_GAP_MAX, sg, gapPeakThr) // 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()) { if (spans.isEmpty()) {
return Result( return Result(
@@ -125,19 +281,37 @@ object DetectLumenFirst {
val lowMean = lowSum / (e - s + 1) val lowMean = lowSum / (e - s + 1)
// py2 method_b — peak_min must satisfy BOTH (a) low_mean + margin, // 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). // (b) >= low_echo_amp (so that wall peaks aren't picked from below T).
val peakMin = maxOf(lowMean + Params.MIN_PEAK_MARGIN, T) val peakMin = maxOf(lowMean + params.minPeakMargin, T)
val ant = WallSelect.selectWallByProminence( var ant = WallSelect.selectWallByProminence(
sg, edge = s, searchWin = Params.PEAK_SEARCH_WIN, sg, edge = s, searchWin = params.peakSearchWin,
peakMin = peakMin, side = WallSelect.Side.ANT, otherEdge = e, peakMin = peakMin, side = WallSelect.Side.ANT, otherEdge = e,
edgeDistDecay = Params.EDGE_DIST_DECAY, maxCandidates = params.maxPeakCandidatesAnt,
valleyStopRise = Params.VALLEY_STOP_RISE, 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( var post = WallSelect.selectWallByProminence(
sg, edge = e, searchWin = Params.PEAK_SEARCH_WIN, sg, edge = e, searchWin = params.peakSearchWin,
peakMin = peakMin, side = WallSelect.Side.POST, otherEdge = s, peakMin = peakMin, side = WallSelect.Side.POST, otherEdge = s,
edgeDistDecay = Params.EDGE_DIST_DECAY, maxCandidates = params.maxPeakCandidatesPost,
valleyStopRise = Params.VALLEY_STOP_RISE, edgeDistDecay = params.edgeDistDecay,
valleyStopRise = params.valleyStopRise,
) )
if (ant == null || post == null) { if (ant == null || post == null) {
return Result( return Result(
@@ -148,15 +322,16 @@ object DetectLumenFirst {
) )
} }
if (post > Params.POST_MAX_IDX) { if (post > params.postMaxIdx) {
val backHalf = ((s + e) / 2) val backHalf = ((s + e) / 2)
val post2 = WallSelect.selectWallByProminence( val post2 = WallSelect.selectWallByProminence(
sg, edge = backHalf, sg, edge = backHalf,
searchWin = Params.POST_MAX_IDX - backHalf, searchWin = params.postMaxIdx - backHalf,
peakMin = peakMin, side = WallSelect.Side.POST, peakMin = peakMin, side = WallSelect.Side.POST,
otherEdge = null, otherEdge = null,
edgeDistDecay = Params.EDGE_DIST_DECAY, maxCandidates = params.maxPeakCandidatesPost,
valleyStopRise = Params.VALLEY_STOP_RISE, edgeDistDecay = params.edgeDistDecay,
valleyStopRise = params.valleyStopRise,
) )
if (post2 == null) { if (post2 == null) {
return Result( return Result(
@@ -170,7 +345,7 @@ object DetectLumenFirst {
} }
val urineLen = post - ant - 1 val urineLen = post - ant - 1
if (urineLen < Params.MIN_URINE_LEN) { if (urineLen < params.minUrineLen) {
return Result( return Result(
mode, raw, sg, cfarThr, T, lowMask, rawSpans, spans, mode, raw, sg, cfarThr, T, lowMask, rawSpans, spans,
ant = null, post = null, lowStart = null, lowEnd = null, ant = null, post = null, lowStart = null, lowEnd = null,
@@ -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.example.medilightv2android.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.example.medilightv2android.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.example.medilightv2android.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)
}
@@ -72,4 +72,55 @@ object Otsu {
} }
return centers[bestT] 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()
}
} }
@@ -29,9 +29,15 @@ object SpanUtils {
} }
/** /**
* Merge spans whose gap ≤ maxGap. Optional peak-guard: * Stage 1 merge — width-bounded amplitude-aware merge. Closes gaps
* if `sg` is provided, gaps where any sample exceeds gapPeakThr are NOT merged * of width ≤ maxGap when the gap-peak amplitude is below gapPeakThr
* (preserves spans separated by a strong peak). * (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( fun mergeCloseSpans(
spans: List<Span>, spans: List<Span>,
@@ -64,4 +70,55 @@ object SpanUtils {
} }
return merged 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,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.example.medilightv2android.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
}
}
@@ -16,8 +16,33 @@ import kotlin.math.min
object WallSelect { object WallSelect {
const val PEAK_SEARCH_WIN = 20 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 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 } enum class Side { ANT, POST }
/** /**
@@ -13,11 +13,26 @@ package com.example.medilightv2android.walldetect.dto
data class BvDispatchResult( data class BvDispatchResult(
val bvMl: Float?, // null = no estimate val bvMl: Float?, // null = no estimate
val rMm: Float?, // equivalent sphere radius (sphere or chord-derived) val rMm: Float?, // equivalent sphere radius (sphere or chord-derived)
val method: String, // FrustumLR | FrustumNoLR | SphereLM | ConeFallback | Verathon | None val method: String, // FrustumLR | FrustumNoLR | ChordMedian | Verathon | None
val confidence: Float, // 0..1 val confidence: Float, // 0..1
val nCenter: Int, // gated center channels (CH0..CH3) val nCenter: Int, // gated center channels (CH0..CH3) AFTER consensus filter
val nLateral: Int, // gated lateral channels (CH4, CH5) val nLateral: Int, // gated lateral channels (CH4, CH5) AFTER consensus filter
val lrRatio: Float, // applied LR/AP ratio (1.0 if no lateral) val lrRatio: Float, // applied LR/AP ratio (1.0 if no lateral)
val warnings: List<String> = emptyList(), val warnings: List<String> = emptyList(),
val sphereCrossCheckBvMl: Float? = null, 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)