水淹模拟以及实现过程（转）

最新推荐文章于 2020-12-23 16:33:06 发布

odexin

最新推荐文章于 2020-12-23 16:33:06 发布

阅读量626

点赞数

分类专栏： skyline

原文链接：https://blog.youkuaiyun.com/sinat_36564005/article/details/109528723

版权

skyline 专栏收录该内容

1 篇文章

订阅专栏

转自：https://blog.youkuaiyun.com/sinat_36564005/article/details/109528723

背景:

在使用skyline开发的过程中,有些功能并不满足实际需求,所以我们就需要自己重写设计以及做一些计算,比如模拟水淹过程,水淹洼地分析,以及水淹面积和体积计算等等.

需求整个计算过程:

水淹算法本身并不难,只是种子填充算法的扩展. 通俗点讲,就是像素之间的连通域问题; 难的是在从鼠标选择区域->切割改区域的地形数据->将三维地形数据降维处理->应用算法处理之后->反映射原始数据->

三维数据面片化->构建棱柱体->计算共面非共面的不规则体积以及面积; 最后才能得到自己想要的功能以及实现,哦,对了,这个结果最后还要可视化;

大家看这个过程,算法本身只是占了一小部分.剩下的用程序展现这个东西才是重点所在,大部分难题都是工程问题.以及开发过程中的种种异常以及流畅体验(异步多线程以及多进程等问题)

先展示一个结果,然后具体聊实现过程

这个结果我和其它主流的分析软件计算的结果对比过,相差不大,只不过其它的软件在计算水淹高度的过程中,可能设置了固定的采样间隔,我没有设置,

所以,可能计算的结果更逼近一些.(虽然这个意义感觉并不大而已)

实现:

数据获取我就不写了,可以用GDAL获取一个tif图像,转三维数据;或者直接从地球上裁剪;具体视情况而异

1.拿到数据之后,我们要将三维数据栅格化.就是

public bool ElevationData2BitmP(ref List<ElevationPoint3D> ListFaces, int width, int height, double waterHeight, out Bitmap a)
{
a = new Bitmap(width, height);
if (width < 1 && height < 1)
{
return false;
}
System.Drawing.Rectangle rect = new System.Drawing.Rectangle(0, 0, a.Width, a.Height);
System.Drawing.Imaging.BitmapData bmpData = a.LockBits(rect, System.Drawing.Imaging.ImageLockMode.ReadWrite,
System.Drawing.Imaging.PixelFormat.Format24bppRgb); //填充数据.
int stride = bmpData.Stride;
int Listnumber = ListFaces.Count;
//查询的步长.
int querySride = width * height / Listnumber;
double lowWaterHeightNumber = 0;
unsafe
{
byte* pIn = (byte*)bmpData.Scan0.ToPointer();
byte* P;
int R, G, B;
double temp = 0;
Random rm = new Random();
for (int y = 0, k = 0; y < a.Height; y++)
{
for (int x = 0; x < a.Width; x++, k += querySride)
{
// double dscal = 0;
P = pIn;
B = P[0];
G = P[1];
R = P[2];
if (k < Listnumber)
{
temp = ListFaces[k].PxyzConvertImg.Z;
if (temp >= waterHeight) //背景色，以及过滤色
{
//dscal = 0;
P[2] = (byte)(73);//R
P[1] = (byte)(221);//G
P[0] = (byte)(93);//B
}
else
{
//dscal = 1;
P[2] = (byte)(0);//R
P[1] = (byte)(0);//G
P[0] = (byte)(255);//B
lowWaterHeightNumber++;
}
}
//P[2] = (byte)(255 * dscal);//R
//P[1] = (byte)(255 * dscal);//G
//P[0] = (byte)(255 * dscal);//B
pIn += 3;
}
pIn += stride - a.Width * 3;
}
}
a.UnlockBits(bmpData);
//得到当前蓝色点.
return true;
}

这里解释一下: 这一部分主要是将三维数据映射成二维uv数据,方便后续计算

通过这一步,就可以实现其数据的可视化. 这一步完成之后,就可以设置一不同高度,显示不同的计算结果

2.计算完成之后,将三维数据网格化

这样可以分析计算每一个小面片的具体参数; 注意其需要依据水淹高度, 在底层构建一个虚假的面.和这一层构成棱柱体.

类似这样.

3.计算,主要是利用海伦公式快速计算

/// <summary>
/// 依据传入的面,以及高. 求取拟柱体体积.
/// </summary> 点的绘制按顺时针计算.
/// <returns></returns>
public double CalculateFitting(ref QuadPrism quadranglePrism)
{
double cylinderVolume = 0;
//以点出发,计算每个点到水平面的距离.
// double currentPlength
//1/6*(S0 + 4S1 + S2)*H
//注意这里的面积,一般都按照梯形面接计算.
//点1 quadrangleFace.item0, item1 ,item2 , item3.
// 从L-> R-> RB-> LB
// item0, item1 item2 item3
// Distance[0] Distance[1] Distance[2] Distance[3]
Point3 currentP = quadranglePrism._bottomFace.pxyz0;
Point3 RightcurrentP = quadranglePrism._bottomFace.pxyz1;
Point3 RightBottomcurrentP = quadranglePrism._bottomFace.pxyz2;
Point3 LeftBottomcurrentP = quadranglePrism._bottomFace.pxyz3;
//SO面积.
double S0h = ComputorFunction.calculateD2Distance(ref currentP, ref LeftBottomcurrentP);
double S0Area = (quadranglePrism.PointToWaterLevelDistance[0] + quadranglePrism.PointToWaterLevelDistance[3]) / 2.0 * S0h;
//S2面积.
double S2h = ComputorFunction.calculateD2Distance(ref RightcurrentP, ref RightBottomcurrentP);
double S2Area = (quadranglePrism.PointToWaterLevelDistance[1] + quadranglePrism.PointToWaterLevelDistance[2]) / 2.0 * S2h;
//S1面积,按照中值划分,可以知道
double S1_upBorder = 0.5 * (quadranglePrism.PointToWaterLevelDistance[0] + quadranglePrism.PointToWaterLevelDistance[1]);
double S1_BottomBorder = 0.5 * (quadranglePrism.PointToWaterLevelDistance[2] + quadranglePrism.PointToWaterLevelDistance[3]);
double S1h = ComputorFunction.calculateD2Distance(ref RightcurrentP, ref RightBottomcurrentP);
double S1Area = (S1_upBorder + S1_BottomBorder) * S1h * 0.5;
//整个体积的高.
double Sh = ComputorFunction.calculateD3Distance(ref currentP, ref RightcurrentP);
//幸普斯规则.数值积分进行拟柱体体积计算.
cylinderVolume = (S0Area + 4 * S1Area + S2Area) / 6.0 * Sh;
return cylinderVolume;
}