From 3b760a3b819a00e58620276179f2d31441464230 Mon Sep 17 00:00:00 2001 From: Zhuangkh Date: Wed, 4 Sep 2024 17:26:43 +0800 Subject: [PATCH] =?UTF-8?q?=E6=8B=9F=E5=90=88?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Commands/Monitor/MonitorViewModel.cs | 177 +++++++++++++++--- uBIMEarthTools/Core/Kriging.cs | 155 +++++++++++++++ uBIMEarthTools/Model/Borehole.cs | 21 +++ uBIMEarthTools/Utils/Common.cs | 32 ++-- 4 files changed, 345 insertions(+), 40 deletions(-) create mode 100644 uBIMEarthTools/Core/Kriging.cs diff --git a/uBIMEarthTools/Commands/Monitor/MonitorViewModel.cs b/uBIMEarthTools/Commands/Monitor/MonitorViewModel.cs index 040c992..cd3dc61 100644 --- a/uBIMEarthTools/Commands/Monitor/MonitorViewModel.cs +++ b/uBIMEarthTools/Commands/Monitor/MonitorViewModel.cs @@ -1,4 +1,5 @@ -using HelixToolkit.Wpf.SharpDX; +using Autodesk.Revit.DB; +using HelixToolkit.Wpf.SharpDX; using SharpDX; using System; using System.Collections.Generic; @@ -78,6 +79,81 @@ namespace uBIMEarthTools.Commands.Monitor } #endregion + #region 拟合 + boreholeList.ForEach(x => x.Interpolation(2)); + + List tris = Common.Delaunay(boreholeList, 10000000 / 304.8); + var topList = boreholeList.Select(x => + { + var p = x.ValueList.First(); + return new Point2D() { X = p.Point.X, Y = p.Point.Y, Value = p.Point.Z }; + }).ToList(); + var topKriging = new IDW(topList, 2); + + var bottomList = boreholeList.Select(x => + { + var p = x.ValueList.Last(); + return new Point2D() { X = p.Point.X, Y = p.Point.Y, Value = p.Point.Z }; + }).ToList(); + var bottomKriging = new IDW(bottomList, 2); + + var typeValue = new Dictionary(); + var valueType = new Dictionary(); + var types = boreholeList.SelectMany(x => x.ValueList.Select(p => p.Type)).Distinct().ToList(); + for (int i = 0; i < types.Count; i++) + { + typeValue.Add(types[i], i + 1); + valueType.Add(i + 1, types[i]); + } + var boreholeKriging = new IDW(boreholeList.SelectMany(b => b.ValueList.Select(p => new Point3D() { X = p.Point.X, Y = p.Point.Y, Z = p.Point.Z, Value = typeValue[p.Type] })).ToList(), 2); + + List boreholes1 = new List(); + //tris.AsParallel().ForAll(item => + //{ + // Func func = (p => + // { + // var topZ = topKriging.Interpolate(new Point2D() { X = p.X, Y = p.Y }); + // var bottomZ = bottomKriging.Interpolate(new Point2D() { X = p.X, Y = p.Y }); + // var plus = (topZ - bottomZ) / 10; + // var geologyLayers = new List(); + // for (int i = 0; i < 10; i++) + // { + // var p3d = new Point3D() { X = p.X, Y = p.Y, Z = topZ - i * plus }; + // var value = boreholeKriging.Interpolate2(p3d); + // geologyLayers.Add(new GeologyLayer(valueType[(int)Math.Round(value)], new Autodesk.Revit.DB.XYZ(p3d.X, p3d.Y, p3d.Z))); + // } + // return new Borehole("", geologyLayers); + // }); + + // var v1 = boreholeList.ElementAt(item[1]).ValueList[0].Point - boreholeList.ElementAt(item[0]).ValueList[0].Point; + // var v2 = boreholeList.ElementAt(item[2]).ValueList[0].Point - boreholeList.ElementAt(item[0]).ValueList[0].Point; + // var originP = boreholeList.ElementAt(item[0]).ValueList[0].Point; + // var triangle = new List() { boreholeList.ElementAt(item[0]).ValueList[0].Point, + // boreholeList.ElementAt(item[1]).ValueList[0].Point, + // boreholeList.ElementAt(item[2]).ValueList[0].Point }; + // var seg = 10; + // v1 = v1 / seg; + // v2 = v2 / seg; + // for (int i = 0; i < seg; i++) + // { + // for (int j = 0; j < seg; j++) + // { + // if (i == j && i == 0) { continue; } + // var xyz = originP + v1 * i + v2 * j; + // if (IsInPolygon(xyz, triangle)) + // { + // boreholes1.Add(func(xyz)); + // } + // } + // } + //}); + //foreach (int[] item in tris) + //{ + + //} + #endregion + boreholeList.AddRange(boreholes1); + #region 计算地质块 List triIndexList = Common.Delaunay(boreholeList, 10000000 / 304.8); List geologyBlocks = new List(); @@ -112,10 +188,13 @@ namespace uBIMEarthTools.Commands.Monitor GeologyCompute.CreateGeologyBlock(boreholes, geologyBlocks, topToBottom); topToBottom = !topToBottom; } + + // break; } #endregion + #region 显示 var colors = new Dictionary(); Random random = new Random(); @@ -124,54 +203,61 @@ namespace uBIMEarthTools.Commands.Monitor colors.Add(item, PhongMaterials.ToColor(random.Next(0, 255) / 255f, random.Next(0, 255) / 255f, random.Next(0, 255) / 255f)); } - foreach (var x in geologyBlocks) + var origin = new SharpDX.Vector3((float)geologyBlocks.First().FLeftUp.X, (float)geologyBlocks.First().FLeftUp.Z, -(float)geologyBlocks.First().FLeftUp.Y); + foreach (var item in geologyBlocks.GroupBy(g => g.Type)) { + var index = 0; List vector3s = new List(); List indices = new List(); - vector3s.Add(new SharpDX.Vector3((float)x.FLeftUp.X, (float)x.FLeftUp.Z, -(float)x.FLeftUp.Y)); //0 - vector3s.Add(new SharpDX.Vector3((float)x.FLeftBottom.X, (float)x.FLeftBottom.Z, -(float)x.FLeftBottom.Y)); //1 - vector3s.Add(new SharpDX.Vector3((float)x.FRightUp.X, (float)x.FRightUp.Z, -(float)x.FRightUp.Y)); //2 - vector3s.Add(new SharpDX.Vector3((float)x.FRightBottom.X, (float)x.FRightBottom.Z, -(float)x.FRightBottom.Y)); //3 - vector3s.Add(new SharpDX.Vector3((float)x.ELeftUp.X, (float)x.ELeftUp.Z, -(float)x.ELeftUp.Y)); //4 - vector3s.Add(new SharpDX.Vector3((float)x.ELeftBottom.X, (float)x.ELeftBottom.Z, -(float)x.ELeftBottom.Y)); //5 - vector3s.Add(new SharpDX.Vector3((float)x.ERightUp.X, (float)x.ERightUp.Z, -(float)x.ERightUp.Y)); //6 - vector3s.Add(new SharpDX.Vector3((float)x.ERightBottom.X, (float)x.ERightBottom.Z, -(float)x.ERightBottom.Y)); //7 + foreach (var x in item) + { + vector3s.Add(new SharpDX.Vector3((float)x.FLeftUp.X, (float)x.FLeftUp.Z, -(float)x.FLeftUp.Y)); //0 + vector3s.Add(new SharpDX.Vector3((float)x.FLeftBottom.X, (float)x.FLeftBottom.Z, -(float)x.FLeftBottom.Y)); //1 + vector3s.Add(new SharpDX.Vector3((float)x.FRightUp.X, (float)x.FRightUp.Z, -(float)x.FRightUp.Y)); //2 + vector3s.Add(new SharpDX.Vector3((float)x.FRightBottom.X, (float)x.FRightBottom.Z, -(float)x.FRightBottom.Y)); //3 + vector3s.Add(new SharpDX.Vector3((float)x.ELeftUp.X, (float)x.ELeftUp.Z, -(float)x.ELeftUp.Y)); //4 + vector3s.Add(new SharpDX.Vector3((float)x.ELeftBottom.X, (float)x.ELeftBottom.Z, -(float)x.ELeftBottom.Y)); //5 + vector3s.Add(new SharpDX.Vector3((float)x.ERightUp.X, (float)x.ERightUp.Z, -(float)x.ERightUp.Y)); //6 + vector3s.Add(new SharpDX.Vector3((float)x.ERightBottom.X, (float)x.ERightBottom.Z, -(float)x.ERightBottom.Y)); //7 - indices.Add(0); indices.Add(2); indices.Add(6); - indices.Add(0); indices.Add(6); indices.Add(4); - indices.Add(0); indices.Add(3); indices.Add(2); - indices.Add(3); indices.Add(0); indices.Add(1); - indices.Add(6); indices.Add(2); indices.Add(3); - indices.Add(6); indices.Add(3); indices.Add(7); + indices.Add(0 + index); indices.Add(2 + index); indices.Add(6 + index); + indices.Add(0 + index); indices.Add(6 + index); indices.Add(4 + index); + indices.Add(0 + index); indices.Add(3 + index); indices.Add(2 + index); + indices.Add(3 + index); indices.Add(0 + index); indices.Add(1 + index); + indices.Add(6 + index); indices.Add(2 + index); indices.Add(3 + index); + indices.Add(6 + index); indices.Add(3 + index); indices.Add(7 + index); - indices.Add(0); indices.Add(4); indices.Add(1); - indices.Add(1); indices.Add(4); indices.Add(5); - indices.Add(4); indices.Add(6); indices.Add(7); - indices.Add(4); indices.Add(7); indices.Add(5); - indices.Add(5); indices.Add(7); indices.Add(1); - indices.Add(1); indices.Add(7); indices.Add(3); + indices.Add(0 + index); indices.Add(4 + index); indices.Add(1 + index); + indices.Add(1 + index); indices.Add(4 + index); indices.Add(5 + index); + indices.Add(4 + index); indices.Add(6 + index); indices.Add(7 + index); + indices.Add(4 + index); indices.Add(7 + index); indices.Add(5 + index); + indices.Add(5 + index); indices.Add(7 + index); indices.Add(1 + index); + indices.Add(1 + index); indices.Add(7 + index); indices.Add(3 + index); + + index = vector3s.Count; + } MeshGeometry3D me = new MeshGeometry3D(); - me.Positions = new Vector3Collection(vector3s); + me.Positions = new Vector3Collection(vector3s.Select(p => p - origin)); me.Indices = new IntCollection(indices); me.Normals = me.CalculateNormals(); - Models.Add(new MeshGeometryModel3D() { Geometry = me, FillMode = SharpDX.Direct3D11.FillMode.Solid, Material = new PhongMaterial() { - AmbientColor = colors[x.Type], - // SpecularColor = new SharpDX.Color4(0.0225f, 0.0225f, 0.0225f, 1.0f), + AmbientColor = colors[item.Key], + // SpecularColor = new SharpDX.Color4(0.0225f, 0.0225f, 0.0225f, 1.0f), //EmissiveColor = new SharpDX.Color4(0.0f, 0.0f, 0.0f, 1.0f), SpecularShininess = 12.8f, - // DiffuseColor = colors[x.Type], + // DiffuseColor = colors[x.Type], } }); - } + #endregion + InitCamera(); } @@ -202,5 +288,40 @@ namespace uBIMEarthTools.Commands.Monitor else return; } + + public static bool IsInPolygon(XYZ checkPoint, List polygonPoints) + { + bool inSide = false; + int pointCount = polygonPoints.Count; + XYZ p1, p2; + for (int i = 0, j = pointCount - 1; i < pointCount; j = i, i++)//第一个点和最后一个点作为第一条线，之后是第一个点和第二个点作为第二条线，之后是第二个点与第三个点，第三个点与第四个点... + { + p1 = polygonPoints[i]; + p2 = polygonPoints[j]; + if (checkPoint.Y < p2.Y) + {//p2在射线之上 + if (p1.Y <= checkPoint.Y) + {//p1正好在射线中或者射线下方 + if ((checkPoint.Y - p1.Y) * (p2.X - p1.X) > (checkPoint.X - p1.X) * (p2.Y - p1.Y))//斜率判断,在P1和P2之间且在P1P2右侧 + { + //射线与多边形交点为奇数时则在多边形之内，若为偶数个交点时则在多边形之外。 + //由于inside初始值为false，即交点数为零。所以当有第一个交点时，则必为奇数，则在内部，此时为inside=(!inside) + //所以当有第二个交点时，则必为偶数，则在外部，此时为inside=(!inside) + inSide = (!inSide); + } + } + } + else if (checkPoint.Y < p1.Y) + { + //p2正好在射线中或者在射线下方，p1在射线上 + if ((checkPoint.Y - p1.Y) * (p2.X - p1.X) < (checkPoint.X - p1.X) * (p2.Y - p1.Y))//斜率判断,在P1和P2之间且在P1P2右侧 + { + inSide = (!inSide); + } + } + } + return inSide; + } + } } diff --git a/uBIMEarthTools/Core/Kriging.cs b/uBIMEarthTools/Core/Kriging.cs new file mode 100644 index 0000000..5c583a1 --- /dev/null +++ b/uBIMEarthTools/Core/Kriging.cs @@ -0,0 +1,155 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; + +namespace uBIMEarthTools.Core +{ + public class IDW where T : Distance + { + private List dataPoints; + private int power; + + public IDW(List dataPoints, int power) + { + this.dataPoints = dataPoints; + this.power = power; + } + + public double Interpolate(T targetPoint) + { + double weightedSum = 0; + double weightSum = 0; + + foreach (var dataPoint in dataPoints) + { + double distance = targetPoint.DistanceTo(dataPoint); + double weight = 1 / Math.Pow(distance, power); + + weightedSum += weight * dataPoint.Value; + weightSum += weight; + } + + return weightedSum / weightSum; + } + + public double Interpolate2(T targetPoint) + { + Dictionary map = new Dictionary(); + foreach (var dataPoint in dataPoints) + { + double distance = targetPoint.DistanceTo(dataPoint); + double weight = 1 / Math.Pow(distance, power); + map[dataPoint.Value] = map.ContainsKey(dataPoint.Value) ? map[dataPoint.Value] + weight : weight; + } + + return map.Aggregate((x,y)=> x.Value > y.Value ? x : y).Key; + } + } + + internal class Kriging where T : Distance + { + private List dataPoints; + private double range; + private double nugget; + + public Kriging(List dataPoints, double range, double nugget) + { + this.dataPoints = dataPoints; + this.range = range; + this.nugget = nugget; + } + + public double Interpolate(T targetPoint) + { + double[] distances = dataPoints.Select(p => targetPoint.DistanceTo(p)).ToArray(); + //range = distances.Max() - distances.Min(); + // Calculate weights based on the spherical variogram model + double[] weights = distances.Select(d => GaussianVariogram(d)).ToArray(); + + // Calculate the interpolated value + double interpolatedValue = 0; + double totalWeight = 0; + + for (int i = 0; i < dataPoints.Count; i++) + { + interpolatedValue += weights[i] * dataPoints[i].Value; + totalWeight += weights[i]; + } + + return interpolatedValue / totalWeight; + } + + private double SphericalVariogram(double h) + { + if (h == 0) + { + return nugget; + } + else if (h <= range) + { + return nugget + (1.5 * Math.Pow(h / range, 3) - 2.5 * Math.Pow(h / range, 2) + 1); + } + else + { + return nugget + 1; + } + } + + private double GaussianVariogram(double h) + { + if (h == 0) + { + return 0; + } + else + { + return nugget + 1 - Math.Exp(-h * h / (range * range)); + } + } + } + + class Point3D : Distance + { + public double X { get; set; } + public double Y { get; set; } + public double Z { get; set; } + + public override double DistanceTo(Point3D other) + { + double dx = this.X - other.X; + double dy = this.Y - other.Y; + double dz = this.Z - other.Z; + + return Math.Sqrt(dx * dx + dy * dy + dz * dz); + } + } + + class Point2D : Distance + { + public double X { get; set; } + public double Y { get; set; } + + public override double DistanceTo(Point2D other) + { + double dx = this.X - other.X; + double dy = this.Y - other.Y; + + return Math.Sqrt(dx * dx + dy * dy ); + } + } + + interface IDistance + { + double DistanceTo(T other); + } + + public abstract class Distance : IDistance + { + public abstract double DistanceTo(T other); + public double Value { get; set; } + } + + +} diff --git a/uBIMEarthTools/Model/Borehole.cs b/uBIMEarthTools/Model/Borehole.cs index 2e8ab82..014d2e4 100644 --- a/uBIMEarthTools/Model/Borehole.cs +++ b/uBIMEarthTools/Model/Borehole.cs @@ -1,4 +1,5 @@ using System.Collections.Generic; +using System.Linq; namespace uBIMEarthTools.Model { @@ -63,5 +64,25 @@ namespace uBIMEarthTools.Model /// 当前地层类型 /// public string CurrentType { get; set; } + + public void Interpolation(double plus) + { + var order = this.ValueList.OrderByDescending(x => x.Point.Z).ToList(); + List newList = new List(); + for (int i = 0; i < order.Count - 1; i++) + { + newList.Add(order[i]); + var z = order[i].Point.Z; + while (z - plus > order[i + 1].Point.Z) + { + z = z - plus; + GeologyLayer geologyLayer = new GeologyLayer(order[i].Type, new Autodesk.Revit.DB.XYZ(order[i].Point.X, order[i].Point.Y, z)); + newList.Add(geologyLayer); + } + } + + this.ValueList = newList; + EndIndex = ValueList.Count - 1; + } } } diff --git a/uBIMEarthTools/Utils/Common.cs b/uBIMEarthTools/Utils/Common.cs index f74e59b..14e3c0c 100644 --- a/uBIMEarthTools/Utils/Common.cs +++ b/uBIMEarthTools/Utils/Common.cs @@ -78,19 +78,27 @@ namespace uBIMEarthTools.Utils /// public static bool IsCollinear(TriangleNet.Topology.Triangle triangle) { - XYZ p0 = new XYZ(triangle.GetVertex(0).X, triangle.GetVertex(0).Y, 0); - XYZ p1 = new XYZ(triangle.GetVertex(1).X, triangle.GetVertex(1).Y, 0); - XYZ p2 = new XYZ(triangle.GetVertex(2).X, triangle.GetVertex(2).Y, 0); - Line line1 = Line.CreateUnbound(p0, p1 - p0); - Line line2 = Line.CreateUnbound(p1, p2 - p1); - Line line3 = Line.CreateUnbound(p2, p0 - p2); - double d1 = line1.Project(p2).Distance / p0.DistanceTo(p1); - double d2 = line2.Project(p0).Distance / p1.DistanceTo(p2); - double d3 = line3.Project(p1).Distance / p2.DistanceTo(p0); - if (Math.Min(Math.Min(d1, d2), d3) < 0.02) + try + { + XYZ p0 = new XYZ(triangle.GetVertex(0).X, triangle.GetVertex(0).Y, 0); + XYZ p1 = new XYZ(triangle.GetVertex(1).X, triangle.GetVertex(1).Y, 0); + XYZ p2 = new XYZ(triangle.GetVertex(2).X, triangle.GetVertex(2).Y, 0); + Line line1 = Line.CreateUnbound(p0, p1 - p0); + Line line2 = Line.CreateUnbound(p1, p2 - p1); + Line line3 = Line.CreateUnbound(p2, p0 - p2); + double d1 = line1.Project(p2).Distance / p0.DistanceTo(p1); + double d2 = line2.Project(p0).Distance / p1.DistanceTo(p2); + double d3 = line3.Project(p1).Distance / p2.DistanceTo(p0); + if (Math.Min(Math.Min(d1, d2), d3) < 0.02) + return true; + else + return false; + } + catch + { return true; - else - return false; + } + } ///