using System; using System.Collections.Generic; using System.ComponentModel; using System.Linq; using System.Runtime.InteropServices; using AI.Common; using AI.Common.Tools; using AI.Common.Log; namespace Reconstruction3DHelper { /// /// 参照opencv里的freeman链码 /// 由于freeman链码里只有0~7,分别表示8个不同方向 /// 但我们的输入点可能不满足上述关系,因此增加了-1和8,分别表示重叠和间距较大的情况 /// internal enum EnumChainCode { // 参考点和当前点重叠 Overlap = -1, Up = 0, UpRight = 1, Right = 2, DownRight = 3, Down = 4, DownLeft = 5, Left = 6, UpLeft = 7, // 当前点不在参考点8邻域内 Gap = 8, } /// /// 线的类型 /// public enum EnumLineType { EightNeighborhood = 0, FourNeighborhood = 1, } /// /// 轮廓相关帮助类 /// public class ContourHelper { /// /// float的精度 /// public const float FloatPrecision = 0.000001f; #region /// /// 参照opencv里的cv::DistanceTypes定义的 /// public enum EnumDistType { User=-1, L1=1, L2=2, C=3, L12=4, Fair=5, Welsch=6, Huber=7, } #endregion #region dllimport [DllImport(@"PlaqueProcessingCpp.dll", CallingConvention = CallingConvention.Cdecl)] public static extern bool FitLine(IntPtr roiVec, int len, float[] matrix, EnumDistType distType, double param, double reps, double aeps); #endregion #region public /// /// 去除轮廓毛刺 /// /// /// /// public static Point2D[] ContourDeburr(Point2D[] points,int filterRadius=3) { // 将轮廓点重采样,使其符合freeman chain code的要求 var resampledPoints = ContourResample(points); // 重采样后点数过少,直接返回 int resampledPointCount = resampledPoints.Length; if (resampledPointCount <= 2) { return resampledPoints; } // 分别计算当前点和前一个点,后一个点的chain code,以此来判断当前点是否为毛刺点 List burrPointIndexes = new List(); for (int ni = 0; ni < resampledPointCount; ni++) { int indexPre = ni; int indexBack = ni; DecrementInCirculation(ref indexPre, resampledPointCount); IncrementInCirculation(ref indexBack, resampledPointCount); if (IsBurrPoint(resampledPoints[indexPre], resampledPoints[ni], resampledPoints[indexBack])) { burrPointIndexes.Add(ni); } } // 如果没有毛刺点,可以直接返回 if (burrPointIndexes.Count <= 0) { return resampledPoints; } // 有毛刺点,需要去除毛刺 // 考虑到毛刺可能突出不止一个像素,因此从每个毛刺点左右分别找与已有轮廓线重合的点 // 只删除与毛刺相邻的重叠点(避免影响其他轮廓交叠位置处) List indexesToRemove = new List(); foreach (var burrPointIndex in burrPointIndexes) { bool preSearch = true; bool backSearch = true; int indexPre = burrPointIndex; int indexBack = burrPointIndex; int searchedCountPre = 0; int searchedCountBack = 0; while (preSearch || backSearch) { // 向前搜索 if (preSearch) { DecrementInCirculation(ref indexPre, resampledPointCount); searchedCountPre += 1; } // 向后搜索 if (backSearch) { IncrementInCirculation(ref indexBack, resampledPointCount); searchedCountBack += 1; } // 判断indexPre处的点,是否和毛刺点后方点相邻,如果相邻,则继续搜索 bool neighborPre = false; for (int ni = 0; ni <= searchedCountBack; ni++) { int index = burrPointIndex; IncrementInCirculation(ref index, resampledPointCount, ni); var chainCodePre = ChainCodeToReferencePoint(resampledPoints[indexPre], resampledPoints[index]); if (chainCodePre != EnumChainCode.Gap) { neighborPre = true; break; } } if (!neighborPre) { preSearch = false; } // 判断indexBack处的点,是否和毛刺点前方的点相邻,如果相邻,则继续搜索 bool neighborBack = false; for (int ni = 0; ni <= searchedCountPre; ni++) { int index = burrPointIndex; DecrementInCirculation(ref index, resampledPointCount, ni); var chainCodeBack = ChainCodeToReferencePoint(resampledPoints[indexBack], resampledPoints[index]); if (chainCodeBack != EnumChainCode.Gap) { neighborBack = true; break; } } if (!neighborBack) { backSearch = false; } if (searchedCountPre + searchedCountBack +1 >= resampledPointCount) { break; } } // 逐个将当前位置处该删除的点添加到总的indexesToRemove里 int indexStart = burrPointIndex; int indexEnd = burrPointIndex; if (searchedCountPre >= 2) { DecrementInCirculation(ref indexStart, resampledPointCount, searchedCountPre-1); } if (searchedCountBack >= 2) { IncrementInCirculation(ref indexEnd, resampledPointCount, searchedCountBack); } int indexRemove = indexStart; while (true) { if (!indexesToRemove.Contains(indexRemove)) { indexesToRemove.Add(indexRemove); } IncrementInCirculation(ref indexRemove, resampledPointCount); if (indexRemove == indexEnd || indexStart == indexEnd) { break; } } } // 没有找到要移除的点,直接返回 if (indexesToRemove.Count <= 0) { return resampledPoints; } // 将待删除的点的序号,排序,从大到小排 indexesToRemove.Sort(); indexesToRemove.Reverse(); List removedPoints = resampledPoints.ToList(); foreach (var index in indexesToRemove) { removedPoints.RemoveAt(index); } // 删除后的点,可能会有空缺,补齐 var pointsDeburred = ContourResample(removedPoints.ToArray()); // 最后再对去毛刺后的点进行滤波 int pointCount = pointsDeburred.Length; if (pointCount <= 2*filterRadius) { return pointsDeburred; } // 得到要滤波的x和y int len = pointCount + 2 * filterRadius; int idx = pointCount - filterRadius; double[] contourX = new double[len]; double[] contourY = new double[len]; for (int ni = 0; ni < len; ni++) { int index = (idx + ni) % pointCount; contourX[ni] = pointsDeburred[index].X; contourY[ni] = pointsDeburred[index].Y; } // 进行一维滤波 double[] contourXFilt = GaussianBlur(contourX, filterRadius); double[] contourYFilt = GaussianBlur(contourY, filterRadius); // 转为输出 Point2D[] pointReturn = new Point2D[pointCount]; for (int ni = 0; ni < pointCount; ni++) { pointReturn[ni].X = (int)contourXFilt[ni + filterRadius]; pointReturn[ni].Y = (int)contourYFilt[ni + filterRadius]; } return pointReturn; } /// /// 轮廓重采样 /// /// /// /// public static Point2D[] ContourResample(Point2D[] points, EnumLineType lineType = EnumLineType.EightNeighborhood) { int pointNum = points.Length; if (pointNum <= 0) { return Array.Empty(); } // 点数过少,直接复制 if (pointNum <= 2) { Point2D[] result = new Point2D[pointNum]; Array.Copy(result, points, pointNum); return result; } // 输入的轮廓点应该是逐像素相连的(8邻域),如果有点重叠或间距较大,则应先重采样成所需格式 List pointsResampled = new List(); Point2D pointPre, pointNi; // 第0个点的前一个点是最后一个点 pointPre = points[pointNum - 1]; for (int ni = 0; ni < pointNum; ni++) { pointNi = points[ni]; // 计算当前点到前一个点之间的连线上所有的点 var linePoints = PointsOnTheLine(pointPre, pointNi, lineType); // 由于前一个点已添加到pointsResampled里,但PointsOnTheLine的返回值是包含起点和终点的,所以从第1个点开始添加 // 当前点和前一个点重叠时,PointsOnTheLine的返回值长度为1,因此不会重复添加重叠点 for (int nj = 1; nj < linePoints.Length; nj++) { pointsResampled.Add(linePoints[nj]); } pointPre = pointNi; } return pointsResampled.ToArray(); } /// /// 循环自减 /// /// /// /// public static void DecrementInCirculation(ref int index, int circleLen, int delta = 1) { delta = delta % circleLen; index -= delta; if (index < 0) { index += circleLen; } } /// /// 循环自增 /// /// /// /// public static void IncrementInCirculation(ref int index, int circleLen, int delta = 1) { delta = delta % circleLen; index += delta; if (index >= circleLen) { index -= circleLen; } } /// /// 一维高斯滤波 /// /// /// /// public static double[] GaussianBlur(double[] datas, int radius) { double[] kernel = GaussianKernel(radius); int len = datas.Length; double[] datasFilt = new double[len]; Array.Copy(datas, datasFilt, len); int kernelLen = kernel.Length; double[] temp = new double[kernelLen]; for (int ni = radius; ni < len - radius; ni++) { double sum = 0; for (int nj = 0; nj < kernelLen; nj++) { temp[nj] = datas[ni + nj - radius] * kernel[nj]; sum += temp[nj]; } datasFilt[ni] = sum; } return datasFilt; } /// /// 获取两点之间的连线 /// 主要参考了opencv的画线算法(该方法用整数的加减计算即可实现,效率较高,属于改进的Bresenham方法) /// (如果pointStart和pointEnd重叠,返回结果只有一个点) /// /// /// /// /// public static Point2D[] PointsOnTheLine(Point2D pointStart, Point2D pointEnd, EnumLineType lineType = EnumLineType.EightNeighborhood) { if (lineType != EnumLineType.EightNeighborhood && lineType != EnumLineType.FourNeighborhood) { throw new ArgumentException("lineType should be 'EightNeighborhood' or 'FourNeighborhood'in this function."); } // 判断是否需要交换xy或交换起始和结束点 int deltaX = 1; int deltaY = 1; int dx = pointEnd.X - pointStart.X; int dy = pointEnd.Y - pointStart.Y; int x = pointStart.X; int y = pointStart.Y; bool vertStartEnd = dx < 0; if (vertStartEnd) { // 交换开始点和结束点 // 注: 这里用if dx<0 则 dx=-dx 来替代Math.abs dx = -dx; dy = -dy; x = pointEnd.X; y = pointEnd.Y; } if (dy < 0) { // 注: 这里用if dy<0 则 dy=-dy 来替代Math.abs dy = -dy; deltaY = -1; } bool vertXY = dy > dx; if (vertXY) { // 交换x和y int temp = dx; dx = dy; dy = temp; temp = deltaX; deltaX = deltaY; deltaY = temp; } // 得到控制循环的一些参数 int error =0, count = 0; int plusErr =0, plusX =0, plusY =0; int minusErr = 0, minusX = 0, minusY = 0; // 注: 这里用 dy+dy 和 dx+dx 代替 2*dy 和 2*dx, 也是为了提高计算效率 switch (lineType) { case EnumLineType.EightNeighborhood: count = dx + 1; error = dx - (dy + dy); plusErr = dx + dx; plusX = 0; plusY = deltaY; minusErr = -(dy + dy); minusX = deltaX; minusY = 0; break; case EnumLineType.FourNeighborhood: count = dx + dy + 1; error = 0; plusErr = (dx + dx) + (dy + dy); plusX = -deltaX; plusY = deltaY; minusErr = -(dy + dy); minusX = deltaX; minusY = 0; break; } if (vertXY) { // 交换plusX和plusY、minusX和minusY int temp = plusX; plusX = plusY; plusY = temp; temp = minusX; minusX = minusY; minusY = temp; } // 循环,得到线上所有点的坐标 Point2D[] points = new Point2D[count]; for (int ni = 0; ni < count; ni++) { int index = ni; if (vertStartEnd) { index = count - 1 - ni; } points[index].X = x; points[index].Y = y; if (error < 0) { y += plusY; error += plusErr; x += plusX; } error += minusErr; x += minusX; y += minusY; } return points; } /// /// 获取一段弧线上的所有点 /// 如果 endAngle大于startAngle 则逆时针绘制,如果endAngle小于startAngle则顺时针绘制 /// 用多边形逼近法(该方法虽然计算量相对较大,但可以较为方便地生成任意的圆弧) /// /// /// /// 起始角度(角度,0~360) /// 角度范围(角度,-360~360,为正则从起始角度逆时针转,为负则顺时针转) /// public static Point2D[] PointsOnTheArc(Point2D pointCenter, float radius, double startAngle, double angleSpan) { // 设置多边形逼近的步长,目前暂定边长约为5个像素 int polygonSideLen = 5; double endAngle = startAngle + angleSpan; // 判断是否需要交换start和end bool vertStartEnd = angleSpan < 0; if (vertStartEnd) { double temp = startAngle; startAngle = endAngle; endAngle = temp; } // 把角度转换到0~360范围内 while (startAngle < 0) { startAngle += 360; endAngle += 360; } while (endAngle > 360) { endAngle -= 360; startAngle -= 360; } if (endAngle - startAngle > 360) { startAngle = 0; endAngle = 360; } // 迭代次数(根据radius和弧度范围,大致确定迭代次数即可) int count = Math.Max(1, Convert.ToInt32(2 * Math.PI * radius * (endAngle - startAngle) / 360 / polygonSideLen)); double alpha = (endAngle - startAngle) / count; double alphaRad = MathHelper.DegreeToRadian(alpha); double startAngleRad = MathHelper.DegreeToRadian(startAngle); double cos = Math.Cos(alphaRad); double sin = Math.Sin(alphaRad); double xi = radius * Math.Cos(startAngleRad); double yi = radius * Math.Sin(startAngleRad); double x, y; List points = new List(); Point2D pointStart = new Point2D(); // 图像坐标上y轴是向下的,所以y应该在pointCenter.Y的基础上减掉yi, 而x方向,因为图像坐标上x轴也是向右,所以是加 pointStart.X = pointCenter.X + (int)xi; pointStart.Y = pointCenter.Y-(int)yi; // 将起始点加到points里 points.Add(pointStart); Point2D pointEnd = new Point2D(); for (int ni = 0; ni < count; ni++) { x = (xi * cos - yi * sin); y = (xi * sin + yi * cos); pointEnd.X = pointCenter.X +(int)x; pointEnd.Y = pointCenter.Y-(int)y; // 两点之间用直线连接 var linePoints = PointsOnTheLine(pointStart, pointEnd); for (int nj = 1; nj < linePoints.Length; nj++) { points.Add(linePoints[nj]); } // 递推 xi = x; yi = y; pointStart = pointEnd; } // 如果交换了起始点和结束点,则翻转一下 if (vertStartEnd) { points.Reverse(); } return points.ToArray(); } /// /// 获取一段椭圆上所有点 /// 如果 endAngle大于startAngle 则逆时针绘制,如果endAngle小于startAngle则顺时针绘制 /// 用多边形逼近法(该方法虽然计算量相对较大,但可以较为方便地生成任意的圆弧) /// /// /// /// /// /// /// /// public static Point2D[] PointsOnTheEllipse(Point2D pointCenter, float width, float height, double angle, double startAngle, double angleSpan) { float halfWidth = width / 2; float halfHeight = height / 2; // 设置多边形逼近的步长,目前暂定边长约为5个像素 int polygonSideLen = 5; // 将angle调整到合适的范围 while (angle < 0) { angle += 360; } while (angle > 360) { angle -= 360; } double endAngle = startAngle + angleSpan; // 判断是否需要交换start和end bool vertStartEnd = angleSpan < 0; if (vertStartEnd) { double temp = startAngle; startAngle = endAngle; endAngle = temp; } // 把角度转换到0~360范围内 while (startAngle < 0) { startAngle += 360; endAngle += 360; } while (endAngle > 360) { endAngle -= 360; startAngle -= 360; } if (endAngle - startAngle > 360) { startAngle = 0; endAngle = 360; } // 迭代次数(根据radius和弧度范围,大致确定迭代次数即可) int count = Math.Max(1, Convert.ToInt32(2 * Math.PI * Math.Sqrt(Math.Pow(halfWidth, 2) + Math.Pow(halfHeight, 2)) * (endAngle - startAngle) / 360 / polygonSideLen)); double alpha = (endAngle - startAngle) / count; double angleRad = MathHelper.DegreeToRadian(angle); double startAngleRad = MathHelper.DegreeToRadian(startAngle); double cosAngle = Math.Cos(angleRad); double sinAngle = Math.Sin(angleRad); double xi = halfWidth * Math.Cos(startAngleRad); double yi = halfHeight * Math.Sin(startAngleRad); double tempAngle,tempRad; List points = new List(); Point2D pointStart = new Point2D(); // 图像坐标上y轴是向下的,所以y应该在pointCenter.Y的基础上减掉yi, 而x方向,因为图像坐标上x轴也是向右,所以是加 pointStart.X = pointCenter.X+(int)(xi * cosAngle - yi * sinAngle); pointStart.Y = pointCenter.Y-(int)(xi * sinAngle + yi * cosAngle); // 将起始点加到points里 points.Add(pointStart); Point2D pointEnd = new Point2D(); for (int ni = 1; ni <= count; ni++) { tempAngle = startAngle + alpha * ni; tempRad = MathHelper.DegreeToRadian(tempAngle); xi = halfWidth * Math.Cos(tempRad); yi = halfHeight * Math.Sin(tempRad); pointEnd.X = pointCenter.X+(int)(xi * cosAngle - yi * sinAngle); pointEnd.Y = pointCenter.Y-(int)(xi * sinAngle + yi * cosAngle); // 两点之间用直线连接 var linePoints = PointsOnTheLine(pointStart, pointEnd); for (int nj = 1; nj < linePoints.Length; nj++) { points.Add(linePoints[nj]); } pointStart = pointEnd; } // 如果交换了起始点和结束点,则翻转一下 if (vertStartEnd) { points.Reverse(); } return points.ToArray(); } /// /// 计算点到直线的距离 /// public static double GetPointToLineDistance(Point2D lineStart, Point2D lineEnd, Point2D point) { if (lineStart.X == lineEnd.X) { return Math.Abs(point.X - lineStart.X); } double slope = (lineEnd.Y - lineStart.Y) / (lineEnd.X - lineStart.X); double intercept = (lineEnd.X * lineStart.Y - lineStart.X * lineEnd.Y) / (lineEnd.X - lineStart.X); return Math.Abs(slope * point.X - point.Y + intercept) / (Math.Sqrt(slope * slope + 1)); } /// /// 计算点到直线的距离 /// /// /// /// public static double GetPointToLineDistance(double[] lineParam, Point2D point) { if (lineParam[1] == 0) { return Math.Abs(point.X - lineParam[2]) / lineParam[0]; } var A = lineParam[0]; var C = lineParam[2]; var a = point.X; var b = point.Y; var n = (A * a + A * A * b + C) / (A * A + 1); var m = a - A * n + A * b; if (Math.Abs(m - a) < 0.00001) { return 0; } var disX = m - a; var disY = n - b; return Math.Sqrt(disX * disX + disY * disY); } /// /// 计算两点间的直线斜率与截距 /// public static void GetSlopeAndIntercept(Point2D FirstPoint, Point2D SeconsPoint, out double k, out double b) { double x1 = FirstPoint.X; double y1 = FirstPoint.Y; double x2 = SeconsPoint.X; double y2 = SeconsPoint.Y; if (x1 == x2) { k = double.NaN; b = double.NaN; return; } if (y1 == y2) { k = 0.0; b = y1; return; } k = ((y1 - y2) / (x1 - x2)); b = y1 - x1 * k; return; } /// /// 射线法判断点是否在多边形内,即计算射线与多边形各边的交点,如果是奇数,则点在多边形内,否则在多边形外 /// 多边形上的点,注意要按照顺序保存 /// 目标点,判断该点是否在多边形内的 /// public static bool IsPointInPolygonRayMethod(Point2D[] polygonPointList, Point2D objectPoint) { if (polygonPointList.Length == 0) { return false; } int num = 0; // 用当前点(i)和前一个点(j)拉一条直线,计算与objectPoint点发出的射线(指向x轴正方向)的交点的个数 int j = polygonPointList.Length - 1; for (int i = 0; i < polygonPointList.Length; i++) { Point2D start = new Point2D(polygonPointList[i].X, polygonPointList[i].Y); Point2D end = new Point2D(polygonPointList[j].X, polygonPointList[j].Y); // 点与多边形的顶点重合 if (((start.X == objectPoint.X) && (start.Y == objectPoint.Y)) || ((end.X == objectPoint.X) && (end.Y == objectPoint.Y))) { return true; } // 判断线段两端点是否在射线两侧 if (((objectPoint.Y >= start.Y) && (objectPoint.Y < end.Y)) || ((objectPoint.Y >= end.Y) && (objectPoint.Y < start.Y))) { // 线段上与射线 Y 坐标相同的点的 X 坐标 var x = start.X + (objectPoint.Y - start.Y) * (end.X - start.X) / (end.Y - start.Y); // 点在多边形的边上 if (objectPoint.X == x) { return true; } // 射线穿过多边形的边界 if (x > objectPoint.X) { ++num; } } j = i; } // 射线穿过多边形边界的次数num为偶数时点在多边形外,返回false,否则返回true return (num % 2 == 0) ? false : true; } /// /// 拟合直线 /// /// /// 是否移除异常点 /// public static bool FitLines(Point2D[] points,out double[] lineParam, bool removeAbnormal = false) { lineParam = new double[3]; try { // 为了避免直接删除输入数组points里的元素,复制一份出来 var pointsCal = new Point2D[points.Length]; Array.Copy(points, pointsCal,points.Length); // 移除异常数据 if (removeAbnormal) { // 原来颈动脉代码里的意思,应该是反复做三次移除异常数据 var removeTimes = 0; while (removeTimes < 3) { pointsCal = RemoveAbnormalData(pointsCal); removeTimes++; } // 因为这里的Point2D,X,Y里存的分别是Point3D的X和Z,所以这里将Point2D按Y排序,即为将投影点仍旧按照其帧顺序排序 // (因为在RemoveAbnormalData里根据点到直线的距离排序后,投影点的帧顺序会被打乱) pointsCal = pointsCal.OrderBy(x => x.Y).ToArray(); } int pointCount = pointsCal.Length; float[] matrix = new float[4]; GCHandle hObject = GCHandle.Alloc(pointsCal, GCHandleType.Pinned); IntPtr pObject = hObject.AddrOfPinnedObject(); // ToDo: 是否需要将拟合直线的参数开放出来(之前颈动脉代码里也是写死的) var ret = FitLine(pObject, pointCount, matrix, EnumDistType.L2, 0, 0.001, 0.001); float x0 = matrix[0]; float x1 = matrix[1]; float x2 = matrix[2]; float x3 = matrix[3]; if (x0 < FloatPrecision) { lineParam[0] = 1; lineParam[1] = 0; lineParam[2] = x2; } else { double k = x1 / x0; double b = x3 - k * x2; lineParam[0] = k; lineParam[1] = -1; lineParam[2] = b; } return true; } catch (Exception excep) { LogHelper.ErrorLog("ContourHelper.FitLine Failed:"+excep); return false; } } #endregion #region private /// /// 计算当前点相对于参考点的chain code /// /// /// /// private static EnumChainCode ChainCodeToReferencePoint(Point2D pointRef, Point2D pointToAdd) { int deltaX = pointToAdd.X - pointRef.X; int deltaY = pointToAdd.Y - pointRef.Y; if (deltaX == 0 && deltaY == 0) { return EnumChainCode.Overlap; } if (deltaX == -1 && deltaY == 0) { return EnumChainCode.Left; } if (deltaX == 1 && deltaY == 0) { return EnumChainCode.Right; } if (deltaX == -1 && deltaY == -1) { return EnumChainCode.UpLeft; } if (deltaX == 0 && deltaY == -1) { return EnumChainCode.Up; } if (deltaX == 1 && deltaY == -1) { return EnumChainCode.UpRight; } if (deltaX == -1 && deltaY == 1) { return EnumChainCode.DownLeft; } if (deltaX == 0 && deltaY == 1) { return EnumChainCode.Down; } if (deltaX == 1 && deltaY == 1) { return EnumChainCode.DownRight; } return EnumChainCode.Gap; } /// /// 计算当前点是否为毛刺点 /// /// /// /// /// private static bool IsBurrPoint(Point2D pointBefore, Point2D pointCur, Point2D pointAfter) { int chainCodeWithBef = (int)ChainCodeToReferencePoint(pointCur, pointBefore); int chainCodeWithAft = (int)ChainCodeToReferencePoint(pointCur, pointAfter); // 只有符合freeman chain code 的相邻点才能这样判断 if (chainCodeWithBef >= 0 && chainCodeWithBef <= 7 && chainCodeWithAft >= 0 && chainCodeWithAft <= 7) { int delta = Math.Abs(chainCodeWithBef - chainCodeWithAft); if (delta<=1 || delta>=7) { return true; } } return false; } /// /// 高斯核 /// /// /// private static double[] GaussianKernel(int radius) { int count = radius * 2 + 1; // sigma 决定高斯核的形态,经验值,选择让核边缘的权重较小,中心权重较大的sigma值 double sigma = radius / 1.5; double[] kernel = new double[count]; double sum = 0; for (int ni = 0; ni < count; ni++) { kernel[ni] = Math.Exp(-Math.Pow(ni - radius, 2) / (2 * Math.Pow(sigma, 2))) / (sigma * Math.Sqrt(2 * Math.PI)); sum += kernel[ni]; } // 除以sum for (int ni = 0; ni < count; ni++) { kernel[ni] /= sum; } return kernel; } /// /// 移除异常数据 /// /// /// private static Point2D[] RemoveAbnormalData(Point2D[] points) { // 如果点数过少,则直接返回 var count = points.Length; if (count < 10) { return points; } if (!FitLines(points, out var line, false)) { return points; } // 计算每个点到直线的距离 var diffs = new double[count]; double sumDiff = 0; for (int ni = 0; ni < count; ni++) { var diff = GetPointToLineDistance(line, points[ni]); sumDiff += diff; diffs[ni] = diff; } // 计算距离值的平均值和标准差 var averageDiff = sumDiff / count; double sumStd = 0; for (int ni = 0; ni < count; ni++) { var diff = diffs[ni]; sumStd += Math.Pow(diff - averageDiff, 2); } sumStd = Math.Sqrt(sumStd / count); var minDiff = averageDiff - 3 * sumStd; var maxDiff = averageDiff + 3 * sumStd; // 删除距离大于平均值加减3倍的标准差的点 for (var ni = 0; ni < count; ni++) { var diff = diffs[ni]; if (diff < minDiff || diff > maxDiff) { diffs[ni] = double.MinValue; } } // 去掉10%的数据(注:颈动脉原来代码里是这样写的,这里先直接照搬,后续再考虑优化) var countRemove = count / 10; var diffList = new List>(); for (var ni = 0; ni < count; ni++) { diffList.Add(new KeyValuePair(ni, diffs[ni])); } var sortedDiffList = diffList.OrderByDescending(o => o.Value); var result = new List(); for (var ni = countRemove; ni < count; ni++) { // 在上一步中,距离大于平均值加减3倍标准差的点(被赋为了double.MinValue),不管在不在10%以内,都应该被舍弃 if (sortedDiffList.ElementAt(ni).Value >=0) { result.Add(points[sortedDiffList.ElementAt(ni).Key]); } } return result.ToArray(); } #endregion } }