private static void PrecomputeTraits(TriangleMesh mesh)
        {


            /// compute all cotans
            /// store each cotan at the foot of each halfedge
            foreach (var ti in mesh.Faces)
            {

                //area calculation
                var l1 = ti.Halfedge.FromVertex.Traits.Position - ti.Halfedge.ToVertex.Traits.Position;
                var l2 = ti.Halfedge.Next.FromVertex.Traits.Position - ti.Halfedge.Next.ToVertex.Traits.Position;
                ti.Traits.Area = Vector3.Cross(l1, l2).Length()/2;

                /// compute cotans
                foreach (var hi in ti.Halfedges)
                {
                    /// TODO_A2 Task 3:
                    /// Implement (extended) cotangents Laplacian matrix (easy)
                    /// here it is recommended to precompute cotans and areas for all faces/halfedges
                    /// You should call it at the beginning of CreateCotanLaplacian and CreateExtendedCotanLaplacian.

                    Vector3 v1 = hi.ToVertex.Traits.Position - hi.FromVertex.Traits.Position;
                    Vector3 v2 = hi.Previous.FromVertex.Traits.Position - hi.Previous.ToVertex.Traits.Position;
                    //v1.Normalize();
                    //v2.Normalize();
                    double angle = Math.Acos(Vector3.Dot(v1, v2) / v1.Length()*v2.Length());
                    if (angle > Math.PI/2) angle = Math.PI - angle;
                    hi.Traits.Cotan = 1.0 / Math.Tan(angle);
                }
            }
            //compute
            foreach (var hi in mesh.Halfedges)
            {
                if(hi.Face == null || hi.Opposite.Face == null) continue;
                //triangle is non-obtuse
                Vector3 v1 = hi.ToVertex.Traits.Position - hi.FromVertex.Traits.Position;
                Vector3 v2 = hi.Previous.FromVertex.Traits.Position - hi.Previous.ToVertex.Traits.Position;

                if (hi.Traits.Cotan >= 0 && hi.Next.Traits.Cotan >= 0 && hi.Previous.Traits.Cotan >= 0)
                {
                    hi.Traits.VoronoiRegionArea = (v1.Length()*hi.Previous.Traits.Cotan
                                                  + v2.Length()*hi.Next.Traits.Cotan)/8.0;
                }
                else if (hi.Traits.Cotan < 0) //triangle is obtuse, angle at halfedge is obtuse
                {
                    hi.Traits.VoronoiRegionArea = hi.Face.Traits.Area/2.0;
                }
                else //Some other angle is obtuse
                {
                    hi.Traits.VoronoiRegionArea = hi.Face.Traits.Area / 4.0;
                }

                //calculate the are for the opposite face of the halfedge
                Vector3 v3 = hi.Opposite.Next.ToVertex.Traits.Position - hi.FromVertex.Traits.Position;
                //triangle is non-obtuse
                if (hi.Opposite.Traits.Cotan >= 0 && hi.Opposite.Next.Traits.Cotan >= 0 && hi.Opposite.Previous.Traits.Cotan >= 0)
                {
                    hi.Traits.VoronoiRegionArea += (v1.Length() * hi.Opposite.Previous.Traits.Cotan
                                                  + v3.Length() * hi.Opposite.Traits.Cotan) / 8.0;

                }
                else if (hi.Traits.Cotan < 0) //triangle is obtuse, angle at halfedge is obtuse
                {
                    hi.Traits.VoronoiRegionArea += hi.Opposite.Face.Traits.Area / 2.0;
                }
                else //Some other angle is obtuse
                {
                    hi.Traits.VoronoiRegionArea += hi.Opposite.Face.Traits.Area / 4.0;
                }


            }

        }