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

转自：https://blog.csdn.net/sinat_36564005/article/details/109528723

背景:

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

需求整个计算过程:

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

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

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

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

 

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

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

实现:

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

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

1. public bool ElevationData2BitmP(ref List<ElevationPoint3D> ListFaces, int width, int height, double waterHeight, out Bitmap a)

2. {

3. a = new Bitmap(width, height);

4. if (width < 1 && height < 1)

5. {

6. return false;

7. }

8.  

9. System.Drawing.Rectangle rect = new System.Drawing.Rectangle(0, 0, a.Width, a.Height);

10. System.Drawing.Imaging.BitmapData bmpData = a.LockBits(rect, System.Drawing.Imaging.ImageLockMode.ReadWrite,

11. System.Drawing.Imaging.PixelFormat.Format24bppRgb); //填充数据.

12. int stride = bmpData.Stride;

13.  

14. int Listnumber = ListFaces.Count;

15. //查询的步长.

16. int querySride = width * height / Listnumber;

17.  

18. double lowWaterHeightNumber = 0;

19. unsafe

20. {

21. byte* pIn = (byte*)bmpData.Scan0.ToPointer();

22. byte* P;

23. int R, G, B;

24. double temp = 0;

25. Random rm = new Random();

26. for (int y = 0, k = 0; y < a.Height; y++)

27. {

28. for (int x = 0; x < a.Width; x++, k += querySride)

29. {

30. // double dscal = 0;

31. P = pIn;

32. B = P[0];

33. G = P[1];

34. R = P[2];

35.  

36. if (k < Listnumber)

37. {

38. temp = ListFaces[k].PxyzConvertImg.Z;

39. if (temp >= waterHeight) //背景色，以及过滤色

40. {

41.  

42. //dscal = 0;

43. P[2] = (byte)(73);//R

44. P[1] = (byte)(221);//G

45. P[0] = (byte)(93);//B

46. }

47. else

48. {

49. //dscal = 1;

50. P[2] = (byte)(0);//R

51. P[1] = (byte)(0);//G

52. P[0] = (byte)(255);//B

53. lowWaterHeightNumber++;

54. }

55. }

56.  

57. //P[2] = (byte)(255 * dscal);//R

58. //P[1] = (byte)(255 * dscal);//G

59. //P[0] = (byte)(255 * dscal);//B

60. pIn += 3;

61. }

62.  

63. pIn += stride - a.Width * 3;

64. }

65.  

66. }

67.  

68. a.UnlockBits(bmpData);

69. //得到当前蓝色点.

70. return true;

71. }

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

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

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

 

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

类似这样.

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

1. /// <summary>

2. /// 依据传入的面,以及高. 求取拟柱体体积.

3. /// </summary> 点的绘制按顺时针计算.

4. /// <returns></returns>

5. public double CalculateFitting(ref QuadPrism quadranglePrism)

6. {

7. double cylinderVolume = 0;

8.  

9. //以点出发,计算每个点到水平面的距离.

10. // double currentPlength

11.  

12. //1/6*(S0 + 4S1 + S2)*H

13. //注意这里的面积,一般都按照梯形面接计算.

14. //点1 quadrangleFace.item0, item1 ,item2 , item3.

15.  

16. // 从L-> R-> RB-> LB

17. // item0, item1 item2 item3

18. // Distance[0] Distance[1] Distance[2] Distance[3]

19. Point3 currentP = quadranglePrism._bottomFace.pxyz0;

20. Point3 RightcurrentP = quadranglePrism._bottomFace.pxyz1;

21. Point3 RightBottomcurrentP = quadranglePrism._bottomFace.pxyz2;

22. Point3 LeftBottomcurrentP = quadranglePrism._bottomFace.pxyz3;

23.  

24. //SO面积.

25. double S0h = ComputorFunction.calculateD2Distance(ref currentP, ref LeftBottomcurrentP);

26. double S0Area = (quadranglePrism.PointToWaterLevelDistance[0] + quadranglePrism.PointToWaterLevelDistance[3]) / 2.0 * S0h;

27.  

28. //S2面积.

29. double S2h = ComputorFunction.calculateD2Distance(ref RightcurrentP, ref RightBottomcurrentP);

30. double S2Area = (quadranglePrism.PointToWaterLevelDistance[1] + quadranglePrism.PointToWaterLevelDistance[2]) / 2.0 * S2h;

31.  

32. //S1面积,按照中值划分,可以知道

33. double S1_upBorder = 0.5 * (quadranglePrism.PointToWaterLevelDistance[0] + quadranglePrism.PointToWaterLevelDistance[1]);

34. double S1_BottomBorder = 0.5 * (quadranglePrism.PointToWaterLevelDistance[2] + quadranglePrism.PointToWaterLevelDistance[3]);

35. double S1h = ComputorFunction.calculateD2Distance(ref RightcurrentP, ref RightBottomcurrentP);

36. double S1Area = (S1_upBorder + S1_BottomBorder) * S1h * 0.5;

37.  

38. //整个体积的高.

39. double Sh = ComputorFunction.calculateD3Distance(ref currentP, ref RightcurrentP);

40. //幸普斯规则.数值积分进行拟柱体体积计算.

41. cylinderVolume = (S0Area + 4 * S1Area + S2Area) / 6.0 * Sh;

42.  

43. return cylinderVolume;

44. }

4.最后使用循环体,累加计算就可以了

 

5.这里总结一下:写的过程中,一定要细致,否则特变容易出一些问题