计算体积

 /// <summary>
        /// 计算体积
        /// </summary>
        /// <param name="TVolumn">填的体积</param>
        /// <param name="WVolumn">挖的体积</param>
        /// <param name="line">分界线</param>
        /// /// <param name="gridPoints">格网角点（含d值）</param>
        public void CalVolumn(out double TVolumn, out double WVolumn,List<PointD[]> gridPoints,out List<PointD> line)
        {
            TVolumn = 0; WVolumn = 0;
            line = new List<PointD>();
            //List<PointD[]> gridPoints = Cald(CalGridPoints(knownPoints, boundary,knownPoint), leveElevation) ;
            
            double area = gridLength * gridWidth;//一个格网的面积
            for (int i = 0; i < gridPoints.Count; i++)
            {
                //填挖与不填不挖的个数
                List<int> TJudge = new List<int>() ;
                List<int> WJudge = new List<int> ();
                List<int> judge = new List<int> ();
                
                for (int j = 0; j < gridPoints[i].Length; j++)
                {
                    if (gridPoints[i][j].d>0)
                    {
                        WJudge.Add(j);
                    }
                    else if (gridPoints[i][j].d<0)
                    {
                        TJudge.Add(j);
                    }
                    else
                    {
                        judge.Add(j);
                    }
                }
                if (judge.Count==4)
                {
                    continue;
                }
                else
                {
                    List<int> adjAngle; int diagonal;
                    //全挖 V=(ha+hb+hc+hd)/4*Sab
                    if (TJudge.Count==0)
                    {
                        WVolumn += (gridPoints[i][0].d + gridPoints[i][1].d + gridPoints[i][2].d + gridPoints[i][3].d) * area / 4.0;
                    }
                    //全填
                    else if (TJudge.Count==4)
                    {
                        TVolumn += (gridPoints[i][0].d + gridPoints[i][1].d + gridPoints[i][2].d + gridPoints[i][3].d) * area / (-4.0);
                    }
                    //填一个角
                    else if(TJudge.Count==1)
                    {
                        //获得邻角与对角
                        FindAdjAngleAndDiagonal(TJudge[0],out adjAngle,out diagonal);
                        //计算体积
                        WVolumn += (2 * (Math.Abs(gridPoints[i][adjAngle[0]].d) + Math.Abs(gridPoints[i][adjAngle[1]].d)) + Math.Abs(gridPoints[i][diagonal].d) - Math.Abs(gridPoints[i][TJudge[0]].d)) / 6 * area;
                        TVolumn += Math.Pow(Math.Abs(gridPoints[i][TJudge[0]].d), 3) / (6 * (Math.Abs(gridPoints[i][TJudge[0]].d) + Math.Abs(gridPoints[i][adjAngle[0]].d)) * (Math.Abs(gridPoints[i][TJudge[0]].d) + Math.Abs(gridPoints[i][adjAngle[1]].d))) * area;
                        //添加分界线的点
                        line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][TJudge[0]], gridPoints[i][adjAngle[0]]));
                        line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][TJudge[0]], gridPoints[i][adjAngle[1]]));
                    }
                    else if (TJudge.Count==2)
                    {
                        //获得邻角与对角
                        FindAdjAngleAndDiagonal(TJudge[0], out adjAngle, out diagonal);
                        //计算体积
                        if (TJudge[1] == diagonal)//对角
                        {
                            WVolumn += (Math.Pow(Math.Abs(gridPoints[i][adjAngle[0]].d), 3) / ((Math.Abs(gridPoints[i][TJudge[0]].d) + Math.Abs(gridPoints[i][adjAngle[0]].d)) * (Math.Abs(gridPoints[i][TJudge[1]].d) + Math.Abs(gridPoints[i][adjAngle[0]].d)))
                                + Math.Pow(Math.Abs(gridPoints[i][adjAngle[1]].d), 3) / ((Math.Abs(gridPoints[i][TJudge[0]].d) + Math.Abs(gridPoints[i][adjAngle[1]].d)) * (Math.Abs(gridPoints[i][TJudge[1]].d) + Math.Abs(gridPoints[i][adjAngle[1]].d)))) * area;
                            TVolumn += (2 * (Math.Abs(gridPoints[i][TJudge[0]].d) + Math.Abs(gridPoints[i][TJudge[1]].d)) - (Math.Abs(gridPoints[i][adjAngle[0]].d) - Math.Abs(gridPoints[i][adjAngle[1]].d))) / 6.0 * area;
                            //添加分界线的点
                            line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][TJudge[0]], gridPoints[i][adjAngle[0]]));
                            line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][TJudge[0]], gridPoints[i][adjAngle[1]]));
                            line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][TJudge[1]], gridPoints[i][adjAngle[0]]));
                            line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][TJudge[1]], gridPoints[i][adjAngle[1]]));
                        }
                        else //邻角
                        {
                            List<int> wAngle = new List<int>();
                            for (int j = 0; j < 4; j++)
                            {
                                if (j != TJudge[0] && j != TJudge[1])
                                {
                                    wAngle.Add(j);
                                }
                            }
                            //计算体积
                            WVolumn += area * Math.Pow(Math.Abs(gridPoints[i][wAngle[0]].d) + Math.Abs(gridPoints[i][wAngle[1]].d), 2) / (4 * (Math.Abs(gridPoints[i][wAngle[0]].d) + Math.Abs(gridPoints[i][wAngle[1]].d) + Math.Abs(gridPoints[i][TJudge[0]].d) + Math.Abs(gridPoints[i][TJudge[1]].d)));
                            TVolumn += area * Math.Pow(Math.Abs(gridPoints[i][TJudge[0]].d) + Math.Abs(gridPoints[i][TJudge[1]].d), 2) / (4 * (Math.Abs(gridPoints[i][wAngle[0]].d) + Math.Abs(gridPoints[i][wAngle[1]].d) + Math.Abs(gridPoints[i][TJudge[0]].d) + Math.Abs(gridPoints[i][TJudge[1]].d)));
                            //添加分界线的点
                            int index = 0;
                            for (int j = 0; j < 2; j++)
                            {
                                if (TJudge[1] != adjAngle[j])
                                {
                                    index = adjAngle[j];
                                }
                            }
                            line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][TJudge[0]], gridPoints[i][index]));
                            line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][TJudge[1]], gridPoints[i][diagonal]));
                        }
                    }
                    else
                    {
                        //获得邻角与对角
                        FindAdjAngleAndDiagonal(WJudge[0], out adjAngle, out diagonal);
                        //计算体积
                        TVolumn += (2 * (Math.Abs(gridPoints[i][adjAngle[0]].d) + Math.Abs(gridPoints[i][adjAngle[1]].d)) + Math.Abs(gridPoints[i][diagonal].d) - Math.Abs(gridPoints[i][WJudge[0]].d)) / 6 * area;
                        WVolumn += Math.Pow(Math.Abs(gridPoints[i][WJudge[0]].d), 3) / (6 * (Math.Abs(gridPoints[i][WJudge[0]].d) + Math.Abs(gridPoints[i][adjAngle[0]].d)) * (Math.Abs(gridPoints[i][WJudge[0]].d) + Math.Abs(gridPoints[i][adjAngle[1]].d))) * area;
                        //添加分界线的点
                        line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][WJudge[0]], gridPoints[i][adjAngle[0]]));
                        line.Add(FindDemarcationPoint(leveElevation, gridPoints[i][WJudge[0]], gridPoints[i][adjAngle[1]]));
                    }
                }
            }
        }