aivfo-tl-3.0/ivf_tl_operate_2.0/control/IvfTl.AutoFocus/Imaging/WellDetector.cs

using System;
using System.Collections.Generic;

namespace IvfTl.AutoFocus.Imaging
{
    /// <summary>well 圆检测结果。</summary>
    public class WellCircle
    {
        public double Cx, Cy;
        public double Radius;
        public double FillRatio;       // 连通域面积 / 外接圆面积，越接近1越圆
        public bool Found;
        public double OffsetX, OffsetY;
        public double OffsetXPct, OffsetYPct;
        public bool Complete;          // well 完整在画面内
        public double AreaPct;         // 连通域占全图比例
        public double Aspect, BoxFill, RConsist; // 调试:宽高比/外接框填充/半径一致性
        public string RejectReason = "";
        public override string ToString()
            => Found
             ? $"圆心({Cx:F0},{Cy:F0}) R={Radius:F0} 偏移X={OffsetXPct:F1}% Y={OffsetYPct:F1}% 完整={Complete} 圆度={FillRatio:F2} 面积={AreaPct:F1}%"
             : "未检出well圆";
    }

    /// <summary>
    /// well 大圆检测（纯C#，无OpenCV）。
    /// well 在背光下是一片中灰亮盘。方法：Otsu阈值 → 最大连通域(BFS) →
    /// 由连通域的质心/面积/外接框算圆心半径，并用圆度+最小尺寸过滤，拒绝噪声小亮斑。
    /// 关键改进：用"最大连通域"而非"全部亮像素质心"，避免分散反光/小亮斑把圆心带偏。
    /// 【移植说明 M2-01】逐字搬自 autofocustool/Imaging/WellDetector.cs，仅改命名空间；
    /// 判据常数(MinAreaPct/aspect/boxFill/rConsist/rFrac 阈值)一字不动(03 §1)。
    /// </summary>
    public static class WellDetector
    {
        /// <summary>检出阈值：连通域至少占全图比例、最小半径(下采样像素)、最低圆度。</summary>
        public static double MinAreaPct = 4.0;
        public static double MinFillRatio = 0.55;

        public static WellCircle Detect(byte[] bgr24, int width, int height, int downsample = 4)
        {
            int dw = width / downsample, dh = height / downsample;
            var gray = new byte[dw * dh];
            int stride = width * 3;
            int[] hist = new int[256];
            for (int y = 0; y < dh; y++)
            {
                int sy = y * downsample;
                for (int x = 0; x < dw; x++)
                {
                    int p = sy * stride + (x * downsample) * 3;
                    byte b = bgr24[p], g = bgr24[p + 1], r = bgr24[p + 2];
                    byte gg = (byte)((b * 29 + g * 150 + r * 77) >> 8);
                    gray[y * dw + x] = gg;
                    hist[gg]++;
                }
            }

            int th = Otsu(hist, dw * dh);
            var res = new WellCircle();

            // 最大连通域(BFS, 4邻域)
            var label = new int[dw * dh];
            int bestCnt = 0, bestX0 = 0, bestY0 = 0, bestX1 = 0, bestY1 = 0;
            long bestSx = 0, bestSy = 0;
            var queue = new Queue<int>();
            for (int idx = 0; idx < dw * dh; idx++)
            {
                if (gray[idx] < th || label[idx] != 0) continue;
                // BFS
                int cnt = 0; long sx = 0, sy = 0;
                int minx = dw, miny = dh, maxx = 0, maxy = 0;
                label[idx] = 1; queue.Enqueue(idx);
                while (queue.Count > 0)
                {
                    int q = queue.Dequeue();
                    int qx = q % dw, qy = q / dw;
                    cnt++; sx += qx; sy += qy;
                    if (qx < minx) minx = qx; if (qx > maxx) maxx = qx;
                    if (qy < miny) miny = qy; if (qy > maxy) maxy = qy;
                    if (qx > 0 && gray[q - 1] >= th && label[q - 1] == 0) { label[q - 1] = 1; queue.Enqueue(q - 1); }
                    if (qx < dw - 1 && gray[q + 1] >= th && label[q + 1] == 0) { label[q + 1] = 1; queue.Enqueue(q + 1); }
                    if (qy > 0 && gray[q - dw] >= th && label[q - dw] == 0) { label[q - dw] = 1; queue.Enqueue(q - dw); }
                    if (qy < dh - 1 && gray[q + dw] >= th && label[q + dw] == 0) { label[q + dw] = 1; queue.Enqueue(q + dw); }
                }
                if (cnt > bestCnt) { bestCnt = cnt; bestSx = sx; bestSy = sy; bestX0 = minx; bestY0 = miny; bestX1 = maxx; bestY1 = maxy; }
            }

            res.AreaPct = 100.0 * bestCnt / (dw * dh);
            if (bestCnt < dw * dh * MinAreaPct / 100.0) { res.Found = false; return res; }

            double cxd = (double)bestSx / bestCnt, cyd = (double)bestSy / bestCnt;
            int bw = bestX1 - bestX0 + 1, bh = bestY1 - bestY0 + 1;

            // 半径：外接框平均半宽 与 面积反推半径
            double rBox = (bw + bh) / 4.0;
            double rArea = Math.Sqrt(bestCnt / Math.PI);
            double rd = (rBox + rArea) / 2;

            // ── 排除斜纹反光带的双判据 ──
            // 1) 外接框宽高比：真well≈1(圆)，反光带细长偏离1
            double aspect = (double)Math.Max(bw, bh) / Math.Max(1, Math.Min(bw, bh));
            // 2) 圆度：连通域面积 / 外接框面积。实心圆≈π/4≈0.785；反光带稀疏远小于此
            double boxFill = (double)bestCnt / (bw * bh);
            // 3) 面积半径与外接框半径一致性：真圆两者接近，反光带差很多
            double rConsist = Math.Min(rArea, rBox) / Math.Max(rArea, rBox);

            res.FillRatio = boxFill; // 用外接框填充率作圆度指标
            res.Aspect = aspect; res.BoxFill = boxFill; res.RConsist = rConsist;

            // 半径合理范围：真well半径约占画面宽15~26%(实测~490/2592=19%)。
            // 斜纹反光带半径常达34%(~870)，用上限拦掉；过小的也排除。
            double rFrac = rd * downsample / width;

            // 真well判据：宽高比<1.5(够方) + 外接框填充>0.6(够实心) + 半径一致>0.65 + 半径占比在[0.10,0.28]
            if (aspect > 1.5 || boxFill < 0.6 || rConsist < 0.65 || rFrac > 0.28 || rFrac < 0.10)
            {
                res.Found = false;
                res.RejectReason = $"aspect={aspect:F2} boxFill={boxFill:F2} rConsist={rConsist:F2} rFrac={rFrac:F2}";
                return res;
            }

            res.Cx = cxd * downsample;
            res.Cy = cyd * downsample;
            res.Radius = rd * downsample;
            res.Found = true;
            res.OffsetX = res.Cx - width / 2.0;
            res.OffsetY = res.Cy - height / 2.0;
            res.OffsetXPct = res.OffsetX / width * 100;
            res.OffsetYPct = res.OffsetY / height * 100;

            // 完整性：外接框是否触画面边(下采样坐标)
            res.Complete = bestX0 > 1 && bestY0 > 1 && bestX1 < dw - 2 && bestY1 < dh - 2;
            return res;
        }

        private static int Otsu(int[] hist, int total)
        {
            long sum = 0;
            for (int i = 0; i < 256; i++) sum += (long)i * hist[i];
            long sumB = 0; int wB = 0; double maxVar = 0; int th = 0;
            for (int t = 0; t < 256; t++)
            {
                wB += hist[t];
                if (wB == 0) continue;
                int wF = total - wB;
                if (wF == 0) break;
                sumB += (long)t * hist[t];
                double mB = (double)sumB / wB;
                double mF = (double)(sum - sumB) / wF;
                double between = (double)wB * wF * (mB - mF) * (mB - mF);
                if (between > maxVar) { maxVar = between; th = t; }
            }
            return th;
        }
    }
}