diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/DetectionSanity.kt b/app/src/main/java/com/example/medilightv2android/walldetect/DetectionSanity.kt new file mode 100644 index 0000000..46a8600 --- /dev/null +++ b/app/src/main/java/com/example/medilightv2android/walldetect/DetectionSanity.kt @@ -0,0 +1,92 @@ +/* + * Copyright (c) 2026 Medithings Co., Ltd. + * Author: Charles KWON OhJun + * 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): List { + val reports = mutableListOf() + 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 + } +} diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/SweepStabilizer.kt b/app/src/main/java/com/example/medilightv2android/walldetect/SweepStabilizer.kt new file mode 100644 index 0000000..1696e4b --- /dev/null +++ b/app/src/main/java/com/example/medilightv2android/walldetect/SweepStabilizer.kt @@ -0,0 +1,137 @@ +/* + * Copyright (c) 2026 Medithings Co., Ltd. + * Author: Charles KWON OhJun + * 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>() + + 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>, + val anomalies: List, + ) + + /** 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>): 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() + 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, 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.deepCopy(): Array = + Array(this.size) { this[it].copyOf() } +} diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/V41Detector.kt b/app/src/main/java/com/example/medilightv2android/walldetect/V41Detector.kt index a06905f..c80cb29 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/V41Detector.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/V41Detector.kt @@ -28,7 +28,7 @@ * ▶ OUTPUT · `dto.DetectionResult` * ─────── * algorithm = "v4_1" - * algorithmVersion = "v4.1.0" + * algorithmVersion = "v4.1.1" * processingMs : Double * perChannel[6] : ChannelResult * ├─ 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.BvEstimation 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.ContrastAux import com.example.medilightv2android.walldetect.algo.DetectLumenFirst 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.SphereFit2Step import com.example.medilightv2android.walldetect.algo.SubsampleRefine @@ -132,7 +134,19 @@ class V41Detector( ) : WallDetector { 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 { val t0 = System.nanoTime() @@ -161,10 +175,16 @@ class V41Detector( // 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 (cr.v41Diag?.gated == true) cr.antIdx else null, - post = if (cr.v41Diag?.gated == true) cr.postIdx else null, + 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( @@ -217,8 +237,19 @@ class V41Detector( val wlDiag = WaveletDenoise.diagnose(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) - val det = DetectLumenFirst.detect(rawInt, DetectLumenFirst.Mode.Adaptive) + // 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 @@ -228,9 +259,14 @@ class V41Detector( // 2.c peaks (sg, all candidates) val peaksAll = PeakDetection.findPeaks1D(sg).toList() - // 3.c subsample refine (parabolic, peak kind) - val antRefined = det.ant?.let { SubsampleRefine.refineParabolic(sg, it, SubsampleRefine.Kind.PEAK) } - val postRefined = det.post?.let { SubsampleRefine.refineParabolic(sg, it, SubsampleRefine.Kind.PEAK) } + // 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() } @@ -242,8 +278,17 @@ class V41Detector( val sContrast: Float = (cR?.contrast ?: 0.0).toFloat() val sContrastTier: String = ContrastAux.tierBand(cR?.contrast) - // 4.b B-mode composite score (RAW envelope) - val sR = BModeScore.score(raw, det.ant, det.post) + // 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(), diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/AnatomicalGate.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/AnatomicalGate.kt index ed6dad8..1f9c74c 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/algo/AnatomicalGate.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/AnatomicalGate.kt @@ -35,12 +35,36 @@ object AnatomicalGate { 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, 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, diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/BModeScore.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/BModeScore.kt index 2c193b0..f894c1e 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/algo/BModeScore.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/BModeScore.kt @@ -41,13 +41,52 @@ 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) - val DEFAULT_WEIGHTS = Weights(far = 0.45, dark = 0.25, ant = 0.15, post = 0.15) + 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) - 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( val sFar: Double, @@ -56,10 +95,16 @@ object BModeScore { val gPost: Double, val lumenMean: 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) @@ -131,6 +176,21 @@ object BModeScore { } 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) @@ -138,26 +198,33 @@ object BModeScore { 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.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), - raw = RawValues(sFar, sDark, gAnt, gPost, lumenMean, farMean, outsideMean), - windows = Windows(lLo = ant + 1, lHi = post, fLo = fLo, fHi = fHi) + 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) ) } } diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/BvEstimation.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/BvEstimation.kt index a9c351f..cc4c995 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/algo/BvEstimation.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/BvEstimation.kt @@ -37,12 +37,24 @@ object BvEstimation { private const val LR_NO_DETECTION = 1.0f private const val AREA_K = PI / 4.0 // (π/4) D² - /** Per-channel detection input — only the integer ant/post matter here. */ - data class Detection(val ant: Int?, val post: Int?) + /** 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. - * Returns a `BvDispatchResult`; never null (always has a method, even "None"). + * 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, @@ -51,37 +63,105 @@ object BvEstimation { 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)" } - val centerIdx = (0..3).filter { walls[it].ant != null && walls[it].post != null } - val lateralIdx = (4..5).filter { walls[it].ant != null && walls[it].post != null } + // ── 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 - 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) - // ── Dispatch (2026-04-29 #3 — py2 _bv_core SSOT, lateral=auxiliary) ─ - // 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) + // ── 3) Dispatch ── val primary = when { - nC == 0 -> none(nC, nL, lrRatio, - if (nL == 0) "no detection" else "lateral only — no SI integration") + 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) - else -> frustum(walls, centerIdx, dps, delayMm, - siDeg, sensorZ, nC, nL, lrRatio, - 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") } - 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("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, + ) } // ────────────────────────────────────────────────────────── diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/ChordConsensus.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/ChordConsensus.kt new file mode 100644 index 0000000..b65f265 --- /dev/null +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/ChordConsensus.kt @@ -0,0 +1,183 @@ +/* + * Copyright (c) 2026 Medithings Co., Ltd. + * Author: Charles KWON OhJun + * 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, + val median: Double?, + val mad: Double?, + val leaderCh: Int?, + val rejected: List + ) + + private fun median(arr: List): Double? { + if (arr.isEmpty()) return null + val s = arr.sorted() + return s[s.size / 2] + } + + private fun mad(arr: List, 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, + dps: Double, + trustScore: Double = DEFAULT_TRUST_SCORE, + chordTol: Double = DEFAULT_CHORD_TOL, + kMad: Double = DEFAULT_K_MAD + ): Result { + val passers = mutableListOf() + 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() + val rejected = mutableListOf() + 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, dps: Double, trusted: Set): MedianBV { + val chords = mutableListOf() + 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) + } +} diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/DetectLumenFirst.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/DetectLumenFirst.kt index ed14b41..9bc1425 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/algo/DetectLumenFirst.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/DetectLumenFirst.kt @@ -17,6 +17,11 @@ import com.example.medilightv2android.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 @@ -24,13 +29,109 @@ object DetectLumenFirst { 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 @@ -69,26 +170,46 @@ object DetectLumenFirst { fun detect(rawAdc: IntArray, mode: Mode = Mode.Fixed()): Result { 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 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(sg) + 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.POST_MAX_IDX + 1, sg.size) + val cap = minOf(params.postMaxIdx + 1, sg.size) val view = DoubleArray(cap) { sg[it] } 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) - .filter { (it.end - it.start + 1) >= Params.LOW_MIN_LEN } - val gapPeakThr = T + Params.GAP_PEAK_MARGIN - val spans = SpanUtils.mergeCloseSpans(rawSpans, Params.MERGE_GAP_MAX, sg, gapPeakThr) + .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( @@ -125,19 +281,37 @@ object DetectLumenFirst { 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.MIN_PEAK_MARGIN, T) + val peakMin = maxOf(lowMean + params.minPeakMargin, T) - val ant = WallSelect.selectWallByProminence( - sg, edge = s, searchWin = Params.PEAK_SEARCH_WIN, + var ant = WallSelect.selectWallByProminence( + sg, edge = s, searchWin = params.peakSearchWin, peakMin = peakMin, side = WallSelect.Side.ANT, otherEdge = e, - edgeDistDecay = Params.EDGE_DIST_DECAY, - valleyStopRise = Params.VALLEY_STOP_RISE, + 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.PEAK_SEARCH_WIN, + sg, edge = e, searchWin = params.peakSearchWin, peakMin = peakMin, side = WallSelect.Side.POST, otherEdge = s, - edgeDistDecay = Params.EDGE_DIST_DECAY, - valleyStopRise = Params.VALLEY_STOP_RISE, + maxCandidates = params.maxPeakCandidatesPost, + edgeDistDecay = params.edgeDistDecay, + valleyStopRise = params.valleyStopRise, ) if (ant == null || post == null) { return Result( @@ -148,15 +322,16 @@ object DetectLumenFirst { ) } - if (post > Params.POST_MAX_IDX) { + if (post > params.postMaxIdx) { val backHalf = ((s + e) / 2) val post2 = WallSelect.selectWallByProminence( sg, edge = backHalf, - searchWin = Params.POST_MAX_IDX - backHalf, + searchWin = params.postMaxIdx - backHalf, peakMin = peakMin, side = WallSelect.Side.POST, otherEdge = null, - edgeDistDecay = Params.EDGE_DIST_DECAY, - valleyStopRise = Params.VALLEY_STOP_RISE, + maxCandidates = params.maxPeakCandidatesPost, + edgeDistDecay = params.edgeDistDecay, + valleyStopRise = params.valleyStopRise, ) if (post2 == null) { return Result( @@ -170,7 +345,7 @@ object DetectLumenFirst { } val urineLen = post - ant - 1 - if (urineLen < Params.MIN_URINE_LEN) { + if (urineLen < params.minUrineLen) { return Result( mode, raw, sg, cfarThr, T, lowMask, rawSpans, spans, ant = null, post = null, lowStart = null, lowEnd = null, diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/ImpulseReject.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/ImpulseReject.kt new file mode 100644 index 0000000..9823b78 --- /dev/null +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/ImpulseReject.kt @@ -0,0 +1,172 @@ +/* + * Copyright (c) 2026 Medithings Co., Ltd. + * Author: Charles KWON OhJun + * 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) + + private fun localMedian(x: DoubleArray, lo: Int, hi: Int): Pair { + 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() + 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() + 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, + val flagged: List + ) + + /** + * 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) + } +} diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/MedianFilter.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/MedianFilter.kt new file mode 100644 index 0000000..55cce9b --- /dev/null +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/MedianFilter.kt @@ -0,0 +1,121 @@ +/* + * Copyright (c) 2026 Medithings Co., Ltd. + * Author: Charles KWON OhJun + * 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 + } +} diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/MorphClose.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/MorphClose.kt new file mode 100644 index 0000000..9b460c4 --- /dev/null +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/MorphClose.kt @@ -0,0 +1,66 @@ +/* + * Copyright (c) 2026 Medithings Co., Ltd. + * Author: Charles KWON OhJun + * 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) +} diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/Otsu.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/Otsu.kt index a76376c..7a6a157 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/algo/Otsu.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/Otsu.kt @@ -72,4 +72,55 @@ object Otsu { } 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() + } } diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/SpanUtils.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/SpanUtils.kt index ac19ad5..8be86db 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/algo/SpanUtils.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/SpanUtils.kt @@ -29,9 +29,15 @@ object SpanUtils { } /** - * Merge spans whose gap ≤ maxGap. Optional peak-guard: - * if `sg` is provided, gaps where any sample exceeds gapPeakThr are NOT merged - * (preserves spans separated by a strong peak). + * 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, @@ -64,4 +70,55 @@ object SpanUtils { } 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, + sg: DoubleArray, + otsuThr: Double, + postMaxIdx: Int, + gapPeakHi: Double? = null + ): List { + 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 + } } diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/StaLta.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/StaLta.kt new file mode 100644 index 0000000..9641c2f --- /dev/null +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/StaLta.kt @@ -0,0 +1,88 @@ +/* + * Copyright (c) 2026 Medithings Co., Ltd. + * Author: Charles KWON OhJun + * 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 + } +} diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/algo/WallSelect.kt b/app/src/main/java/com/example/medilightv2android/walldetect/algo/WallSelect.kt index 73a048b..d1b512a 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/algo/WallSelect.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/algo/WallSelect.kt @@ -16,8 +16,33 @@ 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 } /** diff --git a/app/src/main/java/com/example/medilightv2android/walldetect/dto/BvDispatchResult.kt b/app/src/main/java/com/example/medilightv2android/walldetect/dto/BvDispatchResult.kt index 6e19891..23a0779 100644 --- a/app/src/main/java/com/example/medilightv2android/walldetect/dto/BvDispatchResult.kt +++ b/app/src/main/java/com/example/medilightv2android/walldetect/dto/BvDispatchResult.kt @@ -13,11 +13,26 @@ package com.example.medilightv2android.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 | SphereLM | ConeFallback | Verathon | None + val method: String, // FrustumLR | FrustumNoLR | ChordMedian | Verathon | None val confidence: Float, // 0..1 - val nCenter: Int, // gated center channels (CH0..CH3) - val nLateral: Int, // gated lateral channels (CH4, CH5) + 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 = 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 = emptyList(), + /** Channels rejected by the consensus filter, with reason string. */ + val rejectedChannels: List = 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)