uBIM/uBIMEarthTools
7
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

拟合

Browse Source
This commit is contained in:
Zhuangkh 2024-09-04 17:26:43 +08:00
parent ba930ad5a0
commit 3b760a3b81
4 changed files with 345 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

157
uBIMEarthTools/Commands/Monitor/MonitorViewModel.cs
View File

@ -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<int[]> 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<Point2D>(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<Point2D>(bottomList, 2);
var typeValue = new Dictionary<string, int>();
var valueType = new Dictionary<int, string>();
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<Point3D>(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<Borehole> boreholes1 = new List<Borehole>();
//tris.AsParallel().ForAll(item =>
//{
// Func<XYZ, Borehole> 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<GeologyLayer>();
// 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<XYZ>() { 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<int[]> triIndexList = Common.Delaunay(boreholeList, 10000000 / 304.8);
List<GeologyBlock> geologyBlocks = new List<GeologyBlock>();
@ -112,10 +188,13 @@ namespace uBIMEarthTools.Commands.Monitor
GeologyCompute.CreateGeologyBlock(boreholes, geologyBlocks, topToBottom);
topToBottom = !topToBottom;
}
// break;
}
#endregion
#region 显示
var colors = new Dictionary<string, Color4>();
Random random = new Random();
@ -124,10 +203,14 @@ 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<SharpDX.Vector3> vector3s = new List<SharpDX.Vector3>();
List<int> indices = new List<int>();
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
@ -137,41 +220,44 @@ namespace uBIMEarthTools.Commands.Monitor
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],
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],
}
});
}
#endregion
InitCamera();
}
@ -202,5 +288,40 @@ namespace uBIMEarthTools.Commands.Monitor
else
return;
}
public static bool IsInPolygon(XYZ checkPoint, List<XYZ> 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;
}
}
}

155
uBIMEarthTools/Core/Kriging.cs Normal file
View File

@ -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<T> where T : Distance<T>
{
private List<T> dataPoints;
private int power;
public IDW(List<T> 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<double,double> map = new Dictionary<double,double>();
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<T> where T : Distance<T>
{
private List<T> dataPoints;
private double range;
private double nugget;
public Kriging(List<T> 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<Point3D>
{
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<Point2D>
{
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<T>
{
double DistanceTo(T other);
}
public abstract class Distance<T> : IDistance<T>
{
public abstract double DistanceTo(T other);
public double Value { get; set; }
}
}

21
uBIMEarthTools/Model/Borehole.cs
View File

@ -1,4 +1,5 @@
using System.Collections.Generic;
using System.Linq;
namespace uBIMEarthTools.Model
{
@ -63,5 +64,25 @@ namespace uBIMEarthTools.Model
/// 当前地层类型
/// </summary>
public string CurrentType { get; set; }
public void Interpolation(double plus)
{
var order = this.ValueList.OrderByDescending(x => x.Point.Z).ToList();
List<GeologyLayer> newList = new List<GeologyLayer>();
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;
}
}
}

8
uBIMEarthTools/Utils/Common.cs
View File

@ -77,6 +77,8 @@ namespace uBIMEarthTools.Utils
/// <param name="triangle"></param>
/// <returns></returns>
public static bool IsCollinear(TriangleNet.Topology.Triangle triangle)
{
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);
@ -92,6 +94,12 @@ namespace uBIMEarthTools.Utils
else
return false;
}
catch
{
return true;
}
}
/// <summary>
/// 三角形最长边长度