diff --git a/XP.Calibration/CenterCalibration.cs b/XP.Calibration/CenterCalibration.cs new file mode 100644 index 0000000..a17ced5 --- /dev/null +++ b/XP.Calibration/CenterCalibration.cs @@ -0,0 +1,126 @@ +using System.Drawing; +using Emgu.CV; +using Emgu.CV.Util; +using Emgu.CV.Structure; +using XP.Calibration.Core; +using XP.Calibration.Models; + +namespace XP.Calibration; + +/// +/// 中心校准: 椭圆拟合 + 透视修正 | Center calibration via ellipse fitting + perspective correction +/// +public class CenterCalibration +{ + /// + /// 从投影序列执行中心校准 + /// + /// 投影帧列表 (CV_32F) + /// 几何参数 + /// 校准结果 + public CenterCalibrationResult Calibrate(List projections, GeoParams geo) + { + // 检测各帧球心 + var centers = new List(); + for (int i = 0; i < projections.Count; i++) + { + if (BallDetector.DetectCenter(projections[i], out var c)) + centers.Add(c); + } + + if (centers.Count <= 10) + throw new InvalidOperationException( + $"有效检测帧数不足: {centers.Count}, 至少需要 10 帧"); + + // 椭圆拟合 + var ellipse = FitCentersEllipse(centers); + + // 放大倍率 + double M = geo.DSD / geo.DSO; + + // 从椭圆参数反算倾斜角和 R + double ratio = ellipse.ShortAxis / ellipse.LongAxis; + ratio = Math.Clamp(ratio, 0.0, 1.0); + double alphaRad = Math.Acos(ratio); + double alphaDeg = alphaRad * 180.0 / Math.PI; + + // 长轴 = 2 * R * M / pixelSize → R = longAxis * pixelSize / (2 * M) + double R = ellipse.LongAxis * geo.PixelSize / (2.0 * M); + + // 透视修正 + double deltaPx = R * R * Math.Sin(2.0 * alphaRad) + / (2.0 * geo.DSO * geo.DSO) + * geo.DSD / geo.PixelSize; + + // 长轴方向角 + var rawEllipse = FitEllipseRaw(centers); + double angleDeg = rawEllipse.Angle; + float w = rawEllipse.Size.Width; + float h = rawEllipse.Size.Height; + if (h > w) angleDeg += 90.0f; + + double thetaDeg = 90.0 - angleDeg; + double thetaRad = thetaDeg * Math.PI / 180.0; + + double deltaU = deltaPx * Math.Cos(thetaRad); + double deltaV = deltaPx * (-Math.Sin(thetaRad)); + + double u0 = ellipse.Center.X - deltaU; + double v0 = ellipse.Center.Y - deltaV; + + return new CenterCalibrationResult + { + Ellipse = ellipse, + AlphaDeg = alphaDeg, + R_mm = R, + DeltaPx = deltaPx, + FocalU = u0, + FocalV = v0, + DetectedCenters = centers + }; + } + + /// + /// 从 RAW 文件执行中心校准 (便捷方法) + /// + public CenterCalibrationResult CalibrateFromRaw( + string rawPath, int width, int height, int count, GeoParams geo) + { + var projections = RawReader.ReadFloat32(rawPath, width, height, count); + try + { + return Calibrate(projections, geo); + } + finally + { + foreach (var p in projections) p.Dispose(); + } + } + + /// + /// 对检测到的球心序列做椭圆拟合, 返回长短轴和角度 + /// + private static EllipseResult FitCentersEllipse(List pts) + { + var rotRect = FitEllipseRaw(pts); + float a = rotRect.Size.Width; + float b = rotRect.Size.Height; + + return new EllipseResult + { + Center = rotRect.Center, + LongAxis = Math.Max(a, b), + ShortAxis = Math.Min(a, b), + Angle = rotRect.Angle + }; + } + + /// + /// 调用 OpenCV fitEllipse 返回原始 RotatedRect + /// + private static RotatedRect FitEllipseRaw(List pts) + { + using var vp = new VectorOfPointF(pts.ToArray()); + return CvInvoke.FitEllipse(vp); + } +} diff --git a/XP.Calibration/Core/BallDetector.cs b/XP.Calibration/Core/BallDetector.cs new file mode 100644 index 0000000..de2df2b --- /dev/null +++ b/XP.Calibration/Core/BallDetector.cs @@ -0,0 +1,158 @@ +using System.Drawing; +using Emgu.CV; +using Emgu.CV.CvEnum; +using Emgu.CV.Structure; +using Emgu.CV.Util; + +namespace XP.Calibration.Core; + +/// +/// 球心检测器 | Ball center detector from projection images +/// +public static class BallDetector +{ + /// + /// 是否启用亚像素质心 | Enable sub-pixel centroid + /// + public static bool EnableSubPixel { get; set; } = true; + + /// + /// 从单帧投影中检测球心 (自适应阈值 + 亚像素质心法) + /// + public static bool DetectCenter(Mat srcFloat, out PointF center) + { + center = new PointF(-1, -1); + + // 归一化到 0-255 + using var img8 = new Mat(); + CvInvoke.Normalize(srcFloat, img8, 0, 255, NormType.MinMax); + img8.ConvertTo(img8, DepthType.Cv8U); + CvInvoke.GaussianBlur(img8, img8, new Size(3, 3), 0.8); + + // 自适应阈值 + using var bin = new Mat(); + CvInvoke.AdaptiveThreshold(img8, bin, 255, + AdaptiveThresholdType.MeanC, ThresholdType.Binary, 31, -5); + + // 查找轮廓 + using var contours = new VectorOfVectorOfPoint(); + CvInvoke.FindContours(bin, contours, null, + RetrType.External, ChainApproxMethod.ChainApproxNone); + + if (contours.Size == 0) + return false; + + // 找最大轮廓 + int maxIdx = 0; + double maxArea = 0; + for (int i = 0; i < contours.Size; i++) + { + double area = CvInvoke.ContourArea(contours[i]); + if (area > maxArea) + { + maxArea = area; + maxIdx = i; + } + } + + // 创建掩膜 + using var mask = new Mat(srcFloat.Size, DepthType.Cv8U, 1); + mask.SetTo(new MCvScalar(0)); + CvInvoke.DrawContours(mask, contours, maxIdx, new MCvScalar(255), -1); + + if (EnableSubPixel) + { + center = SubPixelCentroid(srcFloat, mask); + } + else + { + var moments = CvInvoke.Moments(contours[maxIdx]); + center = new PointF( + (float)(moments.M10 / moments.M00), + (float)(moments.M01 / moments.M00)); + } + + return center.X >= 0 && center.Y >= 0; + } + + /// + /// 从单帧投影中检测球心 (Canny 边缘 + 椭圆拟合法) + /// + public static bool DetectCenterByEllipse(Mat srcFloat, out PointF center) + { + center = new PointF(-1, -1); + + using var img8 = new Mat(); + CvInvoke.Normalize(srcFloat, img8, 0, 255, NormType.MinMax); + img8.ConvertTo(img8, DepthType.Cv8U); + CvInvoke.GaussianBlur(img8, img8, new Size(3, 3), 0.8); + + using var edges = new Mat(); + CvInvoke.Canny(img8, edges, 30, 80); + + using var contours = new VectorOfVectorOfPoint(); + CvInvoke.FindContours(edges, contours, null, + RetrType.External, ChainApproxMethod.ChainApproxNone); + + if (contours.Size == 0) + return false; + + int maxIdx = 0; + double maxArea = 0; + for (int i = 0; i < contours.Size; i++) + { + double area = CvInvoke.ContourArea(contours[i]); + if (area > maxArea) + { + maxArea = area; + maxIdx = i; + } + } + + if (contours[maxIdx].Size < 5) + return false; + + using var points = new VectorOfPointF( + Array.ConvertAll(contours[maxIdx].ToArray(), p => new PointF(p.X, p.Y))); + var ellipse = CvInvoke.FitEllipse(points); + center = ellipse.Center; + + return true; + } + + /// + /// 亚像素质心计算 | Sub-pixel centroid using intensity weighting + /// + private static PointF SubPixelCentroid(Mat imgFloat, Mat mask) + { + double sumI = 0, sumX = 0, sumY = 0; + + var imgData = imgFloat.GetData() as float[,]; + var maskData = mask.GetData() as byte[,]; + + if (imgData == null || maskData == null) + return new PointF(-1, -1); + + int rows = imgFloat.Rows; + int cols = imgFloat.Cols; + + for (int y = 0; y < rows; y++) + { + for (int x = 0; x < cols; x++) + { + if (maskData[y, x] != 0) + { + double intensity = imgData[y, x]; + sumI += intensity; + sumX += x * intensity; + sumY += y * intensity; + } + } + } + + if (sumI == 0) + return new PointF(-1, -1); + + return new PointF((float)(sumX / sumI), (float)(sumY / sumI)); + } +} diff --git a/XP.Calibration/Core/DualBallDetector.cs b/XP.Calibration/Core/DualBallDetector.cs new file mode 100644 index 0000000..8dd5d8c --- /dev/null +++ b/XP.Calibration/Core/DualBallDetector.cs @@ -0,0 +1,161 @@ +using System.Drawing; +using Emgu.CV; +using Emgu.CV.CvEnum; +using Emgu.CV.Structure; +using Emgu.CV.Util; + +namespace XP.Calibration.Core; + +/// +/// 双球心检测结果 +/// +public record DualBallResult +{ + /// 第一个球心 (面积较大的) + public PointF Center1 { get; init; } + /// 第二个球心 (面积较小的) + public PointF Center2 { get; init; } + /// 两球心距离 (像素) + public double DistancePx { get; init; } + /// 两球心距离 (mm), 需提供 pixelSize 才有效 + public double DistanceMm { get; init; } + /// 第一个球的轮廓面积 + public double Area1 { get; init; } + /// 第二个球的轮廓面积 + public double Area2 { get; init; } +} + +/// +/// 双球心检测器: 从单帧图像中检测两个球心并计算距离 +/// +public static class DualBallDetector +{ + /// + /// 检测图像中两个最大球体的中心及距离 + /// + /// 输入图像 (CV_32F) + /// 像素物理尺寸 (mm), 0 表示不计算物理距离 + /// 最小轮廓面积阈值, 过滤噪声 + /// 是否启用亚像素质心 + public static DualBallResult? Detect( + Mat srcFloat, + double pixelSize = 0, + double minArea = 50, + bool enableSubPixel = true) + { + // 归一化 + 预处理 + using var img8 = new Mat(); + CvInvoke.Normalize(srcFloat, img8, 0, 255, NormType.MinMax); + img8.ConvertTo(img8, DepthType.Cv8U); + CvInvoke.GaussianBlur(img8, img8, new Size(3, 3), 0.8); + + // 自适应阈值 + using var bin = new Mat(); + CvInvoke.AdaptiveThreshold(img8, bin, 255, + AdaptiveThresholdType.MeanC, ThresholdType.Binary, 31, -5); + + // 查找轮廓 + using var contours = new VectorOfVectorOfPoint(); + CvInvoke.FindContours(bin, contours, null, + RetrType.External, ChainApproxMethod.ChainApproxNone); + + // 按面积排序, 取最大的两个 + var candidates = new List<(int Index, double Area)>(); + for (int i = 0; i < contours.Size; i++) + { + double area = CvInvoke.ContourArea(contours[i]); + if (area >= minArea) + candidates.Add((i, area)); + } + + if (candidates.Count < 2) + return null; + + candidates.Sort((a, b) => b.Area.CompareTo(a.Area)); + var top2 = candidates.Take(2).ToList(); + + // 计算两个球心 + var centers = new PointF[2]; + var areas = new double[2]; + + for (int k = 0; k < 2; k++) + { + int idx = top2[k].Index; + areas[k] = top2[k].Area; + + if (enableSubPixel) + { + using var mask = new Mat(srcFloat.Size, DepthType.Cv8U, 1); + mask.SetTo(new MCvScalar(0)); + CvInvoke.DrawContours(mask, contours, idx, new MCvScalar(255), -1); + centers[k] = SubPixelCentroid(srcFloat, mask); + } + else + { + var moments = CvInvoke.Moments(contours[idx]); + centers[k] = new PointF( + (float)(moments.M10 / moments.M00), + (float)(moments.M01 / moments.M00)); + } + } + + // 计算距离 + double dx = centers[0].X - centers[1].X; + double dy = centers[0].Y - centers[1].Y; + double distPx = Math.Sqrt(dx * dx + dy * dy); + double distMm = pixelSize > 0 ? distPx * pixelSize : 0; + + return new DualBallResult + { + Center1 = centers[0], + Center2 = centers[1], + DistancePx = distPx, + DistanceMm = distMm, + Area1 = areas[0], + Area2 = areas[1] + }; + } + + /// + /// 从 8bit 图像检测 (便捷重载, 支持 CV_8U 输入) + /// + public static DualBallResult? DetectFrom8U( + Mat src8U, + double pixelSize = 0, + double minArea = 50, + bool enableSubPixel = true) + { + using var floatMat = new Mat(); + src8U.ConvertTo(floatMat, DepthType.Cv32F); + return Detect(floatMat, pixelSize, minArea, enableSubPixel); + } + + private static PointF SubPixelCentroid(Mat imgFloat, Mat mask) + { + double sumI = 0, sumX = 0, sumY = 0; + var imgData = imgFloat.GetData() as float[,]; + var maskData = mask.GetData() as byte[,]; + + if (imgData == null || maskData == null) + return new PointF(-1, -1); + + for (int y = 0; y < imgFloat.Rows; y++) + { + for (int x = 0; x < imgFloat.Cols; x++) + { + if (maskData[y, x] != 0) + { + double intensity = imgData[y, x]; + sumI += intensity; + sumX += x * intensity; + sumY += y * intensity; + } + } + } + + if (sumI == 0) + return new PointF(-1, -1); + + return new PointF((float)(sumX / sumI), (float)(sumY / sumI)); + } +} diff --git a/XP.Calibration/Core/Projection.cs b/XP.Calibration/Core/Projection.cs new file mode 100644 index 0000000..ae14aa4 --- /dev/null +++ b/XP.Calibration/Core/Projection.cs @@ -0,0 +1,139 @@ +using XP.Calibration.Models; + +namespace XP.Calibration.Core; + +/// +/// TIGRE 投影模型 (角点法) | TIGRE projection model (corner-based, Siddon) +/// +public static class Projection +{ + private struct Point3D + { + public double X, Y, Z; + public Point3D(double x, double y, double z) { X = x; Y = y; Z = z; } + } + + /// + /// 计算物体原点在探测器上的投影像素坐标 + /// + /// 扫描角 (rad) + /// 旋转轴倾斜角 (rad) + /// 探测器 dYaw/Rx (rad) + /// 探测器 dPitch/Ry (rad) + /// 探测器 dRoll/Rz (rad) + /// 物体 X 偏移 = -R (mm) + /// 探测器 U 偏移 (mm) + /// 探测器 V 偏移 (mm) + /// 几何参数 + /// 输出: U 像素坐标 + /// 输出: V 像素坐标 + public static void ProjectPoint( + double scanAngle, double ay, + double detX, double detY, double detZ, + double offOrigX, double offDetecU, double offDetecV, + GeoParams gp, + out double uPx, out double vPx) + { + double ODD = gp.DSD - gp.DSO; + var S = new Point3D(gp.DSO, 0, 0); + + // 探测器角点初始化 + double py = gp.PixelSize * (0 - (double)gp.NDetecU / 2 + 0.5); + double pz = gp.PixelSize * ((double)gp.NDetecV / 2 - 0.5 - 0); + double puy = gp.PixelSize * (1 - (double)gp.NDetecU / 2 + 0.5); + double pvz = gp.PixelSize * ((double)gp.NDetecV / 2 - 0.5 - 1); + + var P = new Point3D(0, py, pz); + var Pu0 = new Point3D(0, puy, pz); + var Pv0 = new Point3D(0, py, pvz); + + // Step 1: rollPitchYaw (探测器自身旋转) + RollPitchYaw(detZ, detY, detX, ref P); + RollPitchYaw(detZ, detY, detX, ref Pu0); + RollPitchYaw(detZ, detY, detX, ref Pv0); + + // 平移回探测器位置 + P.X -= ODD; Pu0.X -= ODD; Pv0.X -= ODD; + + // Step 2: offDetecU/V 偏移 + P.Y += offDetecU; P.Z += offDetecV; + Pu0.Y += offDetecU; Pu0.Z += offDetecV; + Pv0.Y += offDetecU; Pv0.Z += offDetecV; + + // Step 3: eulerZYZ 扫描旋转 + EulerZYZ(scanAngle, ay, 0, ref P); + EulerZYZ(scanAngle, ay, 0, ref Pu0); + EulerZYZ(scanAngle, ay, 0, ref Pv0); + EulerZYZ(scanAngle, ay, 0, ref S); + + // Step 4: offOrigin 偏移 + P.X -= offOrigX; Pu0.X -= offOrigX; Pv0.X -= offOrigX; + S.X -= offOrigX; + + // deltaU, deltaV + var deltaU = new Point3D(Pu0.X - P.X, Pu0.Y - P.Y, Pu0.Z - P.Z); + var deltaV = new Point3D(Pv0.X - P.X, Pv0.Y - P.Y, Pv0.Z - P.Z); + + // 射线: S → 原点 + var ray = new Point3D(-S.X, -S.Y, -S.Z); + + // 探测器平面法线 + var normal = new Point3D( + deltaU.Y * deltaV.Z - deltaU.Z * deltaV.Y, + deltaU.Z * deltaV.X - deltaU.X * deltaV.Z, + deltaU.X * deltaV.Y - deltaU.Y * deltaV.X); + + // 射线与探测器平面交点 + double pfsDotN = (P.X - S.X) * normal.X + (P.Y - S.Y) * normal.Y + (P.Z - S.Z) * normal.Z; + double rayDotN = ray.X * normal.X + ray.Y * normal.Y + ray.Z * normal.Z; + + if (Math.Abs(rayDotN) < 1e-15) + { + uPx = vPx = -1; + return; + } + + double t = pfsDotN / rayDotN; + var hit = new Point3D(S.X + t * ray.X, S.Y + t * ray.Y, S.Z + t * ray.Z); + var hitP = new Point3D(hit.X - P.X, hit.Y - P.Y, hit.Z - P.Z); + + // 投影到像素坐标 + double dU2 = deltaU.X * deltaU.X + deltaU.Y * deltaU.Y + deltaU.Z * deltaU.Z; + double dV2 = deltaV.X * deltaV.X + deltaV.Y * deltaV.Y + deltaV.Z * deltaV.Z; + double pixelU = (hitP.X * deltaU.X + hitP.Y * deltaU.Y + hitP.Z * deltaU.Z) / dU2; + double pixelV = (hitP.X * deltaV.X + hitP.Y * deltaV.Y + hitP.Z * deltaV.Z) / dV2; + + uPx = pixelU; + vPx = (gp.NDetecV - 1) - pixelV; // TIGRE v 翻转 + } + + /// + /// rollPitchYaw: Rz(dRoll) * Ry(dPitch) * Rx(dYaw) + /// + private static void RollPitchYaw(double dRoll, double dPitch, double dYaw, ref Point3D p) + { + double cr = Math.Cos(dRoll), sr = Math.Sin(dRoll); + double cp = Math.Cos(dPitch), sp = Math.Sin(dPitch); + double cy = Math.Cos(dYaw), sy = Math.Sin(dYaw); + + double x = p.X, y = p.Y, z = p.Z; + p.X = cr * cp * x + (cr * sp * sy - sr * cy) * y + (cr * sp * cy + sr * sy) * z; + p.Y = sr * cp * x + (sr * sp * sy + cr * cy) * y + (sr * sp * cy - cr * sy) * z; + p.Z = -sp * x + cp * sy * y + cp * cy * z; + } + + /// + /// eulerZYZ: Rz(alpha) * Ry(theta) * Rz(psi) + /// + private static void EulerZYZ(double alpha, double theta, double psi, ref Point3D p) + { + double ca = Math.Cos(alpha), sa = Math.Sin(alpha); + double ct = Math.Cos(theta), st = Math.Sin(theta); + double cp = Math.Cos(psi), sp = Math.Sin(psi); + + double x = p.X, y = p.Y, z = p.Z; + p.X = (ca * ct * cp - sa * sp) * x + (-ca * ct * sp - sa * cp) * y + ca * st * z; + p.Y = (sa * ct * cp + ca * sp) * x + (-sa * ct * sp + ca * cp) * y + sa * st * z; + p.Z = -st * cp * x + st * sp * y + ct * z; + } +} diff --git a/XP.Calibration/Core/RawReader.cs b/XP.Calibration/Core/RawReader.cs new file mode 100644 index 0000000..b9498db --- /dev/null +++ b/XP.Calibration/Core/RawReader.cs @@ -0,0 +1,92 @@ +using System.Drawing; +using Emgu.CV; +using Emgu.CV.CvEnum; + +namespace XP.Calibration.Core; + +/// +/// RAW 投影数据读取器 | RAW projection data reader +/// +public static class RawReader +{ + /// + /// 从 RAW 文件读取 float32 投影序列 + /// + /// 文件路径 + /// 图像宽度 + /// 图像高度 + /// 投影帧数 + /// Mat 列表 (CV_32F) + public static List ReadFloat32(string path, int width, int height, int count) + { + var projections = new List(count); + int frameBytes = width * height * sizeof(float); + + using var stream = new FileStream(path, FileMode.Open, FileAccess.Read); + byte[] buffer = new byte[frameBytes]; + + for (int i = 0; i < count; i++) + { + int bytesRead = stream.Read(buffer, 0, frameBytes); + if (bytesRead < frameBytes) + throw new EndOfStreamException( + $"帧 {i} 数据不足: 期望 {frameBytes} 字节, 实际 {bytesRead} 字节"); + + var mat = new Mat(height, width, DepthType.Cv32F, 1); + System.Runtime.InteropServices.Marshal.Copy(buffer, 0, mat.DataPointer, frameBytes); + projections.Add(mat); + } + + return projections; + } + + /// + /// 读取单张 TIFF 图像, 返回 CV_32F 灰度 Mat + /// + /// TIFF 文件路径 + /// Mat (CV_32F, 单通道) + public static Mat ReadTiff(string path) + { + if (!File.Exists(path)) + throw new FileNotFoundException($"TIFF 文件不存在: {path}"); + + // 以原始深度读取 (支持 8/16/32bit TIFF) + var src = CvInvoke.Imread(path, ImreadModes.AnyDepth | ImreadModes.Grayscale); + if (src.IsEmpty) + throw new InvalidOperationException($"无法读取 TIFF: {path}"); + + // 统一转为 CV_32F + if (src.Depth != DepthType.Cv32F) + { + var dst = new Mat(); + src.ConvertTo(dst, DepthType.Cv32F); + src.Dispose(); + return dst; + } + + return src; + } + + /// + /// 批量读取多张 TIFF 图像 (按文件名排序) + /// + /// TIFF 文件所在目录 + /// 搜索模式, 默认 "*.tif" + /// Mat 列表 (CV_32F) + public static List ReadTiffSequence(string directory, string pattern = "*.tif") + { + if (!Directory.Exists(directory)) + throw new DirectoryNotFoundException($"目录不存在: {directory}"); + + var files = Directory.GetFiles(directory, pattern) + .Concat(Directory.GetFiles(directory, pattern + "f")) // 同时匹配 .tif 和 .tiff + .Distinct() + .OrderBy(f => f) + .ToList(); + + if (files.Count == 0) + throw new FileNotFoundException($"目录中未找到 TIFF 文件: {directory}"); + + return files.Select(ReadTiff).ToList(); + } +} diff --git a/XP.Calibration/FullCalibration.cs b/XP.Calibration/FullCalibration.cs new file mode 100644 index 0000000..0f6fe9c --- /dev/null +++ b/XP.Calibration/FullCalibration.cs @@ -0,0 +1,363 @@ +using System.Drawing; +using Emgu.CV; +using Emgu.CV.CvEnum; +using Emgu.CV.Structure; +using Emgu.CV.Util; +using XP.Calibration.Core; +using XP.Calibration.Models; + +namespace XP.Calibration; + +/// +/// 完整几何校准: TIGRE 投影模型 + LM 优化 +/// Full geometric calibration: TIGRE projection model + Levenberg-Marquardt optimization +/// +public class FullCalibration +{ + // 参数索引常量 + private const int IDX_AY = 0; + private const int IDX_DET_Y = 1; + private const int IDX_DET_X = 2; + private const int IDX_DET_Z = 3; + private const int IDX_R = 4; + private const int BASE_PARAMS = 5; + + // 参数边界 + private static readonly double[] LowerBase = { -Math.PI / 2, -Math.PI / 3, -Math.PI / 3, -Math.PI, 0.1 }; + private static readonly double[] UpperBase = { Math.PI / 2, Math.PI / 3, Math.PI / 3, Math.PI, 100 }; + + /// + /// 进度回调 | Progress callback (iteration, rms, params) + /// + public Action? OnProgress { get; set; } + + /// + /// 从投影序列执行完整几何校准 + /// + public FullCalibrationResult Calibrate( + List projections, GeoParams geo, FullCalibrationOptions? options = null) + { + options ??= new FullCalibrationOptions(); + + // 检测各帧球心 + var thetas = new List(); + var uObs = new List(); + var vObs = new List(); + var pts = new List(); + + for (int i = 0; i < projections.Count; i++) + { + if (BallDetector.DetectCenter(projections[i], out var c)) + { + thetas.Add((double)i / (projections.Count - 1) * 2.0 * Math.PI); + uObs.Add(c.X); + vObs.Add(c.Y); + pts.Add(c); + } + } + + if (pts.Count <= 10) + throw new InvalidOperationException($"有效检测帧数不足: {pts.Count}"); + + // 椭圆拟合估算初值 + using var vp = new VectorOfPointF(pts.ToArray()); + var ell = CvInvoke.FitEllipse(vp); + double longAxis = Math.Max(ell.Size.Width, ell.Size.Height); + double shortAxis = Math.Min(ell.Size.Width, ell.Size.Height); + double mag = geo.DSD / geo.DSO; + double rInit = longAxis * geo.PixelSize / (2.0 * mag); + double ayInit = Math.Acos(Math.Min(1.0, shortAxis / longAxis)); + + // 确定参数个数和索引 + int np = BASE_PARAMS; + int idxOffDU = -1, idxOffDV = -1, idxDSO = -1, idxDSD = -1; + + if (options.OptimizeDetectorOffset) + { + idxOffDU = np; idxOffDV = np + 1; np += 2; + } + if (options.OptimizeDistances) + { + idxDSO = np; idxDSD = np + 1; np += 2; + } + + // 构造初始参数 + var p = new double[np]; + p[IDX_AY] = ayInit; + p[IDX_DET_Y] = 0; + p[IDX_DET_X] = 0; + p[IDX_DET_Z] = 0; + p[IDX_R] = rInit; + if (options.OptimizeDetectorOffset) + { + p[idxOffDU] = 0; p[idxOffDV] = 0; + } + if (options.OptimizeDistances) + { + p[idxDSO] = geo.DSO; p[idxDSD] = geo.DSD; + } + + // LM 优化 + var ctx = new LmContext(thetas, uObs, vObs, geo, np, + options.OptimizeDetectorOffset, options.OptimizeDistances, + idxOffDU, idxOffDV, idxDSO, idxDSD); + + bool converged = SolveLM(p, ctx, options.MaxIterations, options.Tolerance); + + // 计算最终 RMS + var finalRes = ComputeResiduals(p, ctx); + double rms = Math.Sqrt(finalRes.Sum(r => r * r) / finalRes.Length); + + double dsoOut = options.OptimizeDistances ? p[idxDSO] : geo.DSO; + double dsdOut = options.OptimizeDistances ? p[idxDSD] : geo.DSD; + + return new FullCalibrationResult + { + AyDeg = p[IDX_AY] * 180.0 / Math.PI, + DetYDeg = p[IDX_DET_Y] * 180.0 / Math.PI, + DetXDeg = p[IDX_DET_X] * 180.0 / Math.PI, + DetZDeg = p[IDX_DET_Z] * 180.0 / Math.PI, + R_mm = p[IDX_R], + OffDetecU_mm = options.OptimizeDetectorOffset ? p[idxOffDU] : 0, + OffDetecV_mm = options.OptimizeDetectorOffset ? p[idxOffDV] : 0, + DSO_mm = dsoOut, + DSD_mm = dsdOut, + RmsPx = rms, + Converged = converged + }; + } + + /// + /// 从 RAW 文件执行完整几何校准 (便捷方法) + /// + public FullCalibrationResult CalibrateFromRaw( + string rawPath, int width, int height, int count, + GeoParams geo, FullCalibrationOptions? options = null) + { + var projections = RawReader.ReadFloat32(rawPath, width, height, count); + try + { + return Calibrate(projections, geo, options); + } + finally + { + foreach (var p in projections) p.Dispose(); + } + } + + #region LM Solver + + private record LmContext( + List Thetas, List UObs, List VObs, + GeoParams Geo, int NP, + bool OptOffset, bool OptDist, + int IdxOffDU, int IdxOffDV, int IdxDSO, int IdxDSD); + + private static double[] ComputeResiduals(double[] p, LmContext ctx) + { + double ay = p[IDX_AY], detY = p[IDX_DET_Y]; + double detX = p[IDX_DET_X], detZ = p[IDX_DET_Z], R = p[IDX_R]; + double offU = ctx.OptOffset ? p[ctx.IdxOffDU] : 0; + double offV = ctx.OptOffset ? p[ctx.IdxOffDV] : 0; + + var gp = new GeoParams + { + DSO = ctx.OptDist ? p[ctx.IdxDSO] : ctx.Geo.DSO, + DSD = ctx.OptDist ? p[ctx.IdxDSD] : ctx.Geo.DSD, + PixelSize = ctx.Geo.PixelSize, + NDetecU = ctx.Geo.NDetecU, + NDetecV = ctx.Geo.NDetecV + }; + + int N = ctx.Thetas.Count; + var res = new double[2 * N]; + for (int i = 0; i < N; i++) + { + Projection.ProjectPoint( + ctx.Thetas[i], ay, detX, detY, detZ, + -R, offU, offV, gp, + out double u, out double v); + res[2 * i] = u - ctx.UObs[i]; + res[2 * i + 1] = v - ctx.VObs[i]; + } + return res; + } + + private static double[,] ComputeJacobian(double[] p, LmContext ctx) + { + int nObs = 2 * ctx.Thetas.Count; + int np = ctx.NP; + var J = new double[nObs, np]; + double eps = 1e-7; + + for (int j = 0; j < np; j++) + { + double[] pp = (double[])p.Clone(); + double[] pm = (double[])p.Clone(); + double step = Math.Max(eps, Math.Abs(p[j]) * eps); + pp[j] = p[j] + step; + pm[j] = p[j] - step; + + var rp = ComputeResiduals(pp, ctx); + var rm = ComputeResiduals(pm, ctx); + + double denom = 2.0 * step; + for (int i = 0; i < nObs; i++) + J[i, j] = (rp[i] - rm[i]) / denom; + } + return J; + } + + private static void ClampParams(double[] p, LmContext ctx) + { + for (int i = 0; i < BASE_PARAMS; i++) + p[i] = Math.Clamp(p[i], LowerBase[i], UpperBase[i]); + + if (ctx.OptOffset) + { + p[ctx.IdxOffDU] = Math.Clamp(p[ctx.IdxOffDU], -50, 50); + p[ctx.IdxOffDV] = Math.Clamp(p[ctx.IdxOffDV], -50, 50); + } + if (ctx.OptDist) + { + p[ctx.IdxDSO] = Math.Clamp(p[ctx.IdxDSO], 50, 1000); + p[ctx.IdxDSD] = Math.Clamp(p[ctx.IdxDSD], 100, 2000); + if (p[ctx.IdxDSD] <= p[ctx.IdxDSO]) + p[ctx.IdxDSD] = p[ctx.IdxDSO] + 10; + } + } + + private bool SolveLM(double[] p, LmContext ctx, int maxIter, double tol) + { + double lambda = 1e-3; + var res = ComputeResiduals(p, ctx); + int m = res.Length; + int np = ctx.NP; + double cost = res.Sum(r => r * r); + + for (int iter = 0; iter < maxIter; iter++) + { + var J = ComputeJacobian(p, ctx); + + // JtJ = J^T * J, Jtr = J^T * r + var JtJ = new double[np, np]; + var Jtr = new double[np]; + for (int i = 0; i < np; i++) + { + for (int j = i; j < np; j++) + { + double sum = 0; + for (int k = 0; k < m; k++) + sum += J[k, i] * J[k, j]; + JtJ[i, j] = sum; + JtJ[j, i] = sum; + } + double s = 0; + for (int k = 0; k < m; k++) + s += J[k, i] * res[k]; + Jtr[i] = s; + } + + // A = JtJ + lambda * diag(1 + JtJ) + var A = new double[np, np]; + var b = new double[np]; + Array.Copy(JtJ, A, JtJ.Length); + for (int i = 0; i < np; i++) + { + A[i, i] += lambda * (1.0 + JtJ[i, i]); + b[i] = -Jtr[i]; + } + + // 解线性方程组 (Cholesky) + if (!SolveLinear(A, b, np, out var delta)) + { + lambda *= 10; + continue; + } + + // 尝试更新 + var np2 = new double[np]; + for (int i = 0; i < np; i++) + np2[i] = p[i] + delta[i]; + ClampParams(np2, ctx); + + var newRes = ComputeResiduals(np2, ctx); + double newCost = newRes.Sum(r => r * r); + + if (newCost < cost) + { + double improvement = cost - newCost; + Array.Copy(np2, p, np); + res = newRes; + cost = newCost; + lambda *= 0.3; + if (lambda < 1e-12) lambda = 1e-12; + + OnProgress?.Invoke(iter, Math.Sqrt(cost / m), p); + + if (improvement < tol) + return true; + } + else + { + lambda *= 10; + if (lambda > 1e12) lambda = 1e12; + } + } + return false; + } + + /// + /// 简单的 Cholesky 分解求解 Ax = b + /// + private static bool SolveLinear(double[,] A, double[] b, int n, out double[] x) + { + x = new double[n]; + var L = new double[n, n]; + + // Cholesky: A = L * L^T + for (int i = 0; i < n; i++) + { + for (int j = 0; j <= i; j++) + { + double sum = 0; + for (int k = 0; k < j; k++) + sum += L[i, k] * L[j, k]; + + if (i == j) + { + double val = A[i, i] - sum; + if (val <= 0) return false; + L[i, j] = Math.Sqrt(val); + } + else + { + L[i, j] = (A[i, j] - sum) / L[j, j]; + } + } + } + + // 前代: L * y = b + var y = new double[n]; + for (int i = 0; i < n; i++) + { + double sum = 0; + for (int k = 0; k < i; k++) + sum += L[i, k] * y[k]; + y[i] = (b[i] - sum) / L[i, i]; + } + + // 回代: L^T * x = y + for (int i = n - 1; i >= 0; i--) + { + double sum = 0; + for (int k = i + 1; k < n; k++) + sum += L[k, i] * x[k]; + x[i] = (y[i] - sum) / L[i, i]; + } + + return true; + } + + #endregion +} diff --git a/XP.Calibration/Models/CalibrationModels.cs b/XP.Calibration/Models/CalibrationModels.cs new file mode 100644 index 0000000..f0b9360 --- /dev/null +++ b/XP.Calibration/Models/CalibrationModels.cs @@ -0,0 +1,52 @@ +using System.Drawing; + +namespace XP.Calibration.Models; + +/// +/// 椭圆拟合结果 | Ellipse fitting result +/// +public record EllipseResult +{ + public PointF Center { get; init; } + public float LongAxis { get; init; } + public float ShortAxis { get; init; } + public float Angle { get; init; } +} + +/// +/// 几何参数 | Geometry parameters for CT system +/// +public class GeoParams +{ + /// 焦点到旋转中心距离 (mm) | Distance Source to Origin + public double DSO { get; set; } + /// 焦点到探测器距离 (mm) | Distance Source to Detector + public double DSD { get; set; } + /// 探测器像素大小 (mm) | Detector pixel size + public double PixelSize { get; set; } + /// 探测器水平像素数 | Detector horizontal pixel count + public int NDetecU { get; set; } + /// 探测器垂直像素数 | Detector vertical pixel count + public int NDetecV { get; set; } +} + +/// +/// 中心校准结果 | Center calibration result +/// +public record CenterCalibrationResult +{ + /// 椭圆拟合结果 + public EllipseResult Ellipse { get; init; } = null!; + /// 倾斜角 (度) | Tilt angle in degrees + public double AlphaDeg { get; init; } + /// 反算半径 R (mm) + public double R_mm { get; init; } + /// 透视偏移量 (像素) | Perspective offset in pixels + public double DeltaPx { get; init; } + /// 修正后焦点投影 U 坐标 + public double FocalU { get; init; } + /// 修正后焦点投影 V 坐标 + public double FocalV { get; init; } + /// 各帧检测到的球心坐标 + public List DetectedCenters { get; init; } = new(); +} diff --git a/XP.Calibration/Models/FullCalibrationOptions.cs b/XP.Calibration/Models/FullCalibrationOptions.cs new file mode 100644 index 0000000..f78a50d --- /dev/null +++ b/XP.Calibration/Models/FullCalibrationOptions.cs @@ -0,0 +1,16 @@ +namespace XP.Calibration.Models; + +/// +/// 完整几何校准选项 | Options for full geometric calibration +/// +public class FullCalibrationOptions +{ + /// 是否优化探测器偏移 offDetecU/V + public bool OptimizeDetectorOffset { get; set; } = false; + /// 是否优化 DSO/DSD 距离 + public bool OptimizeDistances { get; set; } = false; + /// LM 最大迭代次数 + public int MaxIterations { get; set; } = 5000; + /// 收敛阈值 + public double Tolerance { get; set; } = 1e-16; +} diff --git a/XP.Calibration/Models/FullCalibrationResult.cs b/XP.Calibration/Models/FullCalibrationResult.cs new file mode 100644 index 0000000..a3e7fb6 --- /dev/null +++ b/XP.Calibration/Models/FullCalibrationResult.cs @@ -0,0 +1,30 @@ +namespace XP.Calibration.Models; + +/// +/// 完整几何校准结果 | Full geometric calibration result +/// +public record FullCalibrationResult +{ + /// 旋转轴倾斜角 ay (度) | Rotation axis tilt, angles[:,1] + public double AyDeg { get; init; } + /// 探测器 dPitch (度) | Detector pitch, rotDet[:,1] + public double DetYDeg { get; init; } + /// 探测器 dYaw (度) | Detector yaw, rotDet[:,0] + public double DetXDeg { get; init; } + /// 探测器 dRoll (度) | Detector roll, rotDet[:,2] + public double DetZDeg { get; init; } + /// 物体偏移 R (mm) | Object offset radius + public double R_mm { get; init; } + /// 探测器水平偏移 (mm) | Detector U offset + public double OffDetecU_mm { get; init; } + /// 探测器垂直偏移 (mm) | Detector V offset + public double OffDetecV_mm { get; init; } + /// 优化后 DSO (mm) + public double DSO_mm { get; init; } + /// 优化后 DSD (mm) + public double DSD_mm { get; init; } + /// RMS 残差 (像素) | RMS residual in pixels + public double RmsPx { get; init; } + /// 是否收敛 | Whether optimization converged + public bool Converged { get; init; } +} diff --git a/XP.Calibration/README.md b/XP.Calibration/README.md new file mode 100644 index 0000000..5c37c82 --- /dev/null +++ b/XP.Calibration/README.md @@ -0,0 +1,159 @@ +# XP.Calibration + +平面 CT 系统几何校准库，基于 .NET 8 + Emgu.CV，从 C++ OpenCV 实现移植而来。 + +## 功能 + +提供两种校准方法： + +- **中心校准 (CenterCalibration)** — 从投影序列检测球心轨迹，椭圆拟合后反算倾斜角和焦点投影偏移，适用于快速估算 +- **完整几何校准 (FullCalibration)** — 基于 TIGRE 投影模型（角点法），使用 Levenberg-Marquardt 优化器同时优化 5~9 个几何参数 + +## 项目结构 + +``` +XP.Calibration/ +├── Models/ +│ ├── CalibrationModels.cs # EllipseResult, GeoParams, CenterCalibrationResult +│ ├── FullCalibrationResult.cs # 完整校准输出 (各参数 + RMS) +│ └── FullCalibrationOptions.cs # 优化模式选项 +├── Core/ +│ ├── RawReader.cs # RAW float32 投影数据读取 +│ ├── BallDetector.cs # 球心检测 (自适应阈值+亚像素质心 / Canny+椭圆拟合) +│ └── TigreProjection.cs # TIGRE 投影模型 (rollPitchYaw + eulerZYZ) +├── CenterCalibration.cs # 中心校准 +└── FullCalibration.cs # 完整几何校准 (LM 优化) +``` + +## 校准参数 + +完整几何校准支持三种模式： + +| 模式 | 参数数 | 优化参数 | +|------|--------|----------| +| 基础 | 5 | ay, det_y, det_x, det_z, R | +| +探测器偏移 | 7 | 基础 + offDetecU, offDetecV | +| +距离 | 9 | 基础 + offDetecU, offDetecV + DSO, DSD | + +参数含义： + +| 参数 | 说明 | TIGRE 对应 | +|------|------|-----------| +| ay | 旋转轴倾斜角 | angles[:,1] | +| det_x | 探测器 dYaw (Rx) | rotDetector[:,0] | +| det_y | 探测器 dPitch (Ry) | rotDetector[:,1] | +| det_z | 探测器 dRoll (Rz) | rotDetector[:,2] | +| R | 物体偏移半径 | offOrigin[:,2] = -R | +| offDetecU/V | 探测器平移偏移 | offDetector[:,0]/[:,1] | +| DSO | 焦点到旋转中心距离 | geo.DSO | +| DSD | 焦点到探测器距离 | geo.DSD | + +## 使用示例 + +### 中心校准 + +```csharp +var geo = new GeoParams +{ + DSO = 200, DSD = 500, PixelSize = 0.6, + NDetecU = 512, NDetecV = 512 +}; + +var calib = new CenterCalibration(); +var result = calib.CalibrateFromRaw("projections.raw", 512, 512, 360, geo); + +Console.WriteLine($"倾斜角: {result.AlphaDeg:F2}°"); +Console.WriteLine($"R: {result.R_mm:F2} mm"); +Console.WriteLine($"焦点投影: ({result.FocalU:F2}, {result.FocalV:F2})"); +``` + +### 完整几何校准 (5 参数) + +```csharp +var geo = new GeoParams +{ + DSO = 250, DSD = 450, PixelSize = 0.6, + NDetecU = 512, NDetecV = 512 +}; + +var full = new FullCalibration(); +var result = full.CalibrateFromRaw("projections.raw", 512, 512, 360, geo); + +Console.WriteLine($"ay: {result.AyDeg:F4}°"); +Console.WriteLine($"det_z: {result.DetZDeg:F4}°"); +Console.WriteLine($"R: {result.R_mm:F4} mm"); +Console.WriteLine($"RMS: {result.RmsPx:F6} px"); +``` + +### 完整几何校准 (9 参数) + +```csharp +var options = new FullCalibrationOptions +{ + OptimizeDetectorOffset = true, // +offDetecU/V + OptimizeDistances = true, // +DSO/DSD + MaxIterations = 5000, + Tolerance = 1e-16 +}; + +var full = new FullCalibration(); +// 可选: 监听优化进度 +full.OnProgress = (iter, rms, p) => + Console.WriteLine($"Iter {iter}: RMS={rms:F6} px"); + +var result = full.CalibrateFromRaw("projections.raw", 512, 512, 360, geo, options); +``` + +### 直接传入 Mat 列表 + +如果投影数据已经在内存中（比如从相机采集），可以直接传 `List`： + +```csharp +List projections = ...; // 已有的 CV_32F Mat 列表 +var result = new FullCalibration().Calibrate(projections, geo, options); +``` + +### 切换球心检测方法 + +默认使用自适应阈值 + 亚像素质心法。如需使用 Canny + 椭圆拟合法： + +```csharp +// 关闭亚像素质心，改用矩心 +BallDetector.EnableSubPixel = false; + +// 或直接调用椭圆拟合检测 +BallDetector.DetectCenterByEllipse(mat, out var center); +``` + +## 算法说明 + +### 球心检测流程 + +1. 归一化 float32 → 8bit +2. 高斯模糊 (3×3, σ=0.8) +3. 自适应阈值二值化 (blockSize=31, C=-5) +4. 查找最大轮廓 +5. 亚像素质心 (强度加权) 或矩心 + +### TIGRE 投影模型 + +完全按照 TIGRE `computeDeltas_Siddon` 流程实现： + +1. 初始化探测器角点 (P, Pu0, Pv0) 和射线源 S +2. `rollPitchYaw` — 探测器自身旋转 (det_x, det_y, det_z) +3. 探测器偏移 (offDetecU/V) +4. `eulerZYZ` — 扫描旋转 (scanAngle, ay) +5. 物体偏移 (offOrigX = -R) +6. 射线-平面交点 → 像素坐标 + +### LM 优化器 + +- 数值差分雅可比矩阵 (中心差分, ε=1e-7) +- Cholesky 分解求解法方程 +- 自适应阻尼因子 (成功 ×0.3, 失败 ×10) +- 参数边界约束 (clamp) + +## 依赖 + +- .NET 8.0-windows +- Emgu.CV 4.10.0.5680 diff --git a/XP.Calibration/XP.Calibration.csproj b/XP.Calibration/XP.Calibration.csproj new file mode 100644 index 0000000..1b93f25 --- /dev/null +++ b/XP.Calibration/XP.Calibration.csproj @@ -0,0 +1,16 @@ + + + + net8.0-windows + enable + enable + XP.Calibration + XP.Calibration + + + + + + + + diff --git a/XP.Hardware.Detector/bin/Debug/net8.0-windows7.0/XP.Hardware.Detector.deps.json b/XP.Hardware.Detector/bin/Debug/net8.0-windows7.0/XP.Hardware.Detector.deps.json index 3f740dc..54eff8d 100644 --- a/XP.Hardware.Detector/bin/Debug/net8.0-windows7.0/XP.Hardware.Detector.deps.json +++ b/XP.Hardware.Detector/bin/Debug/net8.0-windows7.0/XP.Hardware.Detector.deps.json @@ -1878,7 +1878,10 @@ "Telerik.UI.for.Wpf.NetCore.Xaml": "2024.1.408" }, "runtime": { - "XP.Common.dll": {} + "XP.Common.dll": { + "assemblyVersion": "1.4.16.1", + "fileVersion": "1.4.16.1" + } }, "resources": { "en-US/XP.Common.resources.dll": { diff --git a/XP.ImageProcessing.RoiControl/Controls/MeasureEventArgs.cs b/XP.ImageProcessing.RoiControl/Controls/MeasureEventArgs.cs new file mode 100644 index 0000000..7278450 --- /dev/null +++ b/XP.ImageProcessing.RoiControl/Controls/MeasureEventArgs.cs @@ -0,0 +1,31 @@ +using System.Windows; + +namespace XP.ImageProcessing.RoiControl.Controls +{ + /// 测量完成事件参数 + public class MeasureCompletedEventArgs : RoutedEventArgs + { + public Point P1 { get; } + public Point P2 { get; } + public double Distance { get; } + public int TotalCount { get; } + public string MeasureType { get; set; } + + public MeasureCompletedEventArgs(RoutedEvent routedEvent, Point p1, Point p2, double distance, int totalCount) + : base(routedEvent) + { + P1 = p1; P2 = p2; Distance = distance; TotalCount = totalCount; + } + } + + /// 测量状态变化事件参数 + public class MeasureStatusEventArgs : RoutedEventArgs + { + public string Message { get; } + + public MeasureStatusEventArgs(RoutedEvent routedEvent, string message) : base(routedEvent) + { + Message = message; + } + } +} diff --git a/XP.ImageProcessing.RoiControl/Controls/PolygonRoiCanvas.xaml b/XP.ImageProcessing.RoiControl/Controls/PolygonRoiCanvas.xaml index 409d5d4..be035cd 100644 --- a/XP.ImageProcessing.RoiControl/Controls/PolygonRoiCanvas.xaml +++ b/XP.ImageProcessing.RoiControl/Controls/PolygonRoiCanvas.xaml @@ -17,35 +17,20 @@ - - - - - - - - -