using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Microsoft.Xna.Framework;
using System.Diagnostics;

namespace DOFerVolumeRenderer
{
    public class Slicer
    {
        private const float EPSILON = 0.00001f;
        private Vector3[] m_iPoints = new Vector3[6];
        private float[] m_psuedoAngles = new float[6];
		private ICamera m_camera = null;


		public Slicer(ICamera cam)
		{
			m_camera = cam;
		}



        /// <summary>
        /// Ray Plane Intersection
        /// </summary>
        /// <param name="rayOrigin"> Origin of the Ray Vector3</param>
        /// <param name="rayDir"> Direction of the Ray vector3</param>
        /// <param name="plane"> The plane to test a Intersection with. XNA Plane struct</param>
        /// <param name="outT"> if greater or equal 0 and smaller or equal 1, intersecion point is inbetween raysegment</param>
        /// <returns></returns>
        private bool RayPlaneIntersection(Vector3 rayOrigin, Vector3 rayDir, Plane  plane, out float outT)
        {
            float vd = 
                plane.Normal.X * rayDir.X + 
                plane.Normal.Y * rayDir.Y +
                plane.Normal.Z * rayDir.Z;

            outT =
                -(plane.Normal.X * rayOrigin.X +
                  plane.Normal.Y * rayOrigin.Y +
                  plane.Normal.Z * rayOrigin.Z + 
                  plane.D) / vd;       
            return true;
        }


        /*! Returns the pseudoangle between the line p1 to (infinity, p1.y) and the 
        line from p1 to p2. The pseudoangle has the property that the ordering of 
        points by true angle anround p1 and ordering of points by pseudoangle are the 
        same The result is in the range [0, 4) (or error -1). */
        /// <summary>
        /// 
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <returns></returns>
        private static float PseudoAngle(float p1X, float p1Y, float p2X, float p2Y) {
           Vector2 delta = new Vector2(p2X - p1X, p2Y - p1Y);
           float result;
   
           if ((delta.X == 0) && (delta.Y == 0)) {
              return -1;
           } else {
              result = delta.Y / (Math.Abs(delta.X) + Math.Abs(delta.Y));
      
              if (delta.X < 0.0) {
                 result = 2.0f - result;
              } else {
                 result = 4.0f + result;
              }
      
           }
   
           return result;
        }

        // TODO: Special case for point Count 3
        // Todo: Optimize using NOT the Center (less triangles)
        private void SortPointsAndTesselate(int pointCount, ref Vector3[] intersectionPoints)
        {
            //Stopwatch sWatch = new Stopwatch();
            //sWatch.Start();

            // calculate Center
            Vector3 center = m_iPoints[0];

            for (int i = 1; i < pointCount; i++)
            {
                center += m_iPoints[i];
            }
            center = center / pointCount;

            // Project on XY and calculate Pseudo Angles
//            float[] psuedoAngles = new float[pointCount];
            for (int i = 0; i < 6; i++)
            {
                if (i < pointCount)
                    m_psuedoAngles[i] = PseudoAngle(m_iPoints[i].X, m_iPoints[i].Y, center.X, center.Y);
                else
                    m_psuedoAngles[i] = 1000.0f; // 1000 is always Higher than values generated from pseudo Angles
            }
            
            // sort unsorted
            Array.Sort(m_psuedoAngles, m_iPoints);


            // Tesselate
            //Vector3[] sortTess = new Vector3[3 * pointCount];
            for (int i = 0; i < 6; i++)
            {
                if (i < pointCount)
                {
                    intersectionPoints[i * 3] = center;
                    intersectionPoints[i * 3 + 1] = m_iPoints[i];
                    intersectionPoints[i * 3 + 2] = m_iPoints[(i + 1) % pointCount];
                }
                else
                {
                    intersectionPoints[i * 3] = Vector3.Zero;
                    intersectionPoints[i * 3 + 1] = Vector3.Zero;
                    intersectionPoints[i * 3 + 2] = Vector3.Zero;
                }
            }

        }

        /// <summary>
        ///         __________ bbMax
        ///        /         /|
        ///       /         / |
        ///      /_________/  |     y    z
        ///      |         |  /      | /
        ///      |         | /       |/
        ///      |_________|/        °--- x
        ///     bbMin 
        /// 
        /// </summary>
        /// <param name="plane"></param>
        /// <param name="bbMin"></param>
        /// <param name="bbMax"></param>
        /// <param name="intersectionPoints"></param>
        /// <param name="pointCount"></param>
        public void planeBoxIntersection(Plane plane, Vector3 bbMin, Vector3 bbMax, ref Vector3[] intersectionPoints, out int pointCount)
        {

            pointCount = 0;
            float t;

            //Vector3.Transform(bbMin, m_camera.ViewMatrix);
            //Vector3.Transform(bbMax, m_camera.ViewMatrix);


            //Vector3[] iPoints = new Vector3[6];

            // X - Axis Edges
            Vector3 dir = new Vector3(bbMax.X - bbMin.X, 0.0f, 0.0f);
           
            Vector3 origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform(origin + dir * t,m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMax.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            // Y - Axis Edges
            dir = new Vector3( 0.0f, bbMax.Y - bbMin.Y, 0.0f);

            origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMax.X, bbMin.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMax.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            // Z - Axis Edges
            dir = new Vector3(0.0f, 0.0f, bbMax.Y - bbMin.Y);

            origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMax.X, bbMin.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMax.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            if (pointCount > 0)
                SortPointsAndTesselate(pointCount, ref intersectionPoints);
        }



        /// <summary>
        ///         __________ bbMax
        ///        /         /|
        ///       /         / |
        ///      /_________/  |     y    z
        ///      |         |  /      | /
        ///      |         | /       |/
        ///      |_________|/        °--- x
        ///     bbMin 
        /// 
        /// </summary>
        /// <param name="plane"></param>
        /// <param name="bbMin"></param>
        /// <param name="bbMax"></param>
        /// <param name="intersectionPoints"></param>
        /// <param name="pointCount"></param>
        public void planeBoxIntersection2(Plane plane, Vector3 bbMin, Vector3 bbMax, ref Vector3[] intersectionPoints, out int pointCount)
        {

            pointCount = 0;
            float t;

            //Vector3.Transform(bbMin, m_camera.ViewMatrix);
            //Vector3.Transform(bbMax, m_camera.ViewMatrix);


            //Vector3[] iPoints = new Vector3[6];

            // X - Axis Edges
            Vector3 dir = new Vector3(bbMax.X - bbMin.X, 0.0f, 0.0f);

            Vector3 origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform(origin + dir * t, m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMax.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            // Y - Axis Edges
            dir = new Vector3(0.0f, bbMax.Y - bbMin.Y, 0.0f);

            origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMax.X, bbMin.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMax.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            // Z - Axis Edges
            dir = new Vector3(0.0f, 0.0f, bbMax.Y - bbMin.Y);

            origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMax.X, bbMin.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMin.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            origin = new Vector3(bbMax.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= EPSILON && t <= 1.0f)
                m_iPoints[pointCount++] = Vector3.Transform((origin + dir * t), m_camera.ViewMatrix);

            if (pointCount > 0)
                SortPointsAndTesselate(pointCount, ref intersectionPoints);
        }

    }
}
/*
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Microsoft.Xna.Framework;

namespace DOFerVolumeRenderer
{
    class Slicer
    {

        /// <summary>
        /// Ray Plane Intersection
        /// </summary>
        /// <param name="rayOrigin"> Origin of the Ray Vector3</param>
        /// <param name="rayDir"> Direction of the Ray vector3</param>
        /// <param name="plane"> The plane to test a Intersection with. XNA Plane struct</param>
        /// <param name="outT"> if greater or equal 0 and smaller or equal 1, intersecion point is inbetween raysegment</param>
        /// <returns></returns>
        private static bool RayPlaneIntersection(Vector3 rayOrigin, Vector3 rayDir, Plane  plane, out float outT)
        {
            float vd = 
                plane.Normal.X * rayDir.X + 
                plane.Normal.Y * rayDir.Y +
                plane.Normal.Z * rayDir.Z;

            outT =
                -(plane.Normal.X * rayOrigin.X +
                  plane.Normal.Y * rayOrigin.Y +
                  plane.Normal.Z * rayOrigin.Z + 
                  plane.D) / vd;       
            return true;
        }


        /// <summary>
        /// 
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <returns></returns>
        private static float PseudoAngle(Vector2 p1, Vector2 p2) {
           Vector2 delta = p2 - p1;
           float result;
   
           if ((delta.X == 0) && (delta.Y == 0)) {
              return -1;
           } else {
              result = delta.Y / (Math.Abs(delta.X) + Math.Abs(delta.Y));
      
              if (delta.X < 0.0) {
                 result = 2.0f - result;
              } else {
                 result = 4.0f + result;
              }
      
           }
   
           return result;
        }

        public static void TestPseudeAngle()
        {
            Vector2[] pS = new Vector2[5];
            pS[0] = new Vector2(0.25f, 0);
            pS[1] = new Vector2(0, 0.75f);
            pS[2] = new Vector2(0.75f, 0f);
            pS[3] = new Vector2(1, 0.75f);
            pS[4] = new Vector2(0.5f, 1f);

            Vector2 center = (pS[0] + pS[1] + pS[2] + pS[3] + pS[4]) / pS.Length;

            float[] pFC = new float[5];
            pFC[0] = PseudoAngle(pS[0], center);
            pFC[1] = PseudoAngle(pS[1], center);
            pFC[2] = PseudoAngle(pS[2], center);
            pFC[3] = PseudoAngle(pS[3], center);
            pFC[4] = PseudoAngle(pS[4], center);

            Array.Sort(pFC, pS);

        }

        // TODO: Special case for point Count 3
        // Todo: Optimize using NOT the Center (less triangles)
        private static Vector3[] SortPointsAndTesselate(Vector3[] unsorted, int pointCount)
        {
            // calculate Center
            Vector3 center = unsorted[0];
            for (int i = 1; i < pointCount; i++)
            {
                center += unsorted[i];
            }
            center = center / pointCount;

            // Project on XY and calculate Pseudo Angles
            float[] psuedoAngles = new float[pointCount];
            for (int i = 0; i < pointCount; i++)
            {
                psuedoAngles[i] = PseudoAngle(new Vector2(unsorted[i].X, unsorted[i].Y), new Vector2(center.X, center.Y));
            }
            
            // sort unsorted
            Array.Sort(psuedoAngles, unsorted);

            // Tesselate
            Vector3[] sortTess = new Vector3[3 * pointCount];
            for (int i = 0; i < pointCount; i++)
            {
                sortTess[i * 3] = center;
                sortTess[i * 3 + 1] = unsorted[i];
                sortTess[i * 3 + 2] = unsorted[(i + 1) % pointCount];
            }

            return sortTess;
        }

        /// <summary>
        ///         __________ bbMax
        ///        /         /|
        ///       /         / |
        ///      /_________/  |     y    z
        ///      |         |  /      | /
        ///      |         | /       |/
        ///      |_________|/        °--- x
        ///     bbMin 
        /// 
        /// </summary>
        /// <param name="plane"></param>
        /// <param name="bbMin"></param>
        /// <param name="bbMax"></param>
        /// <param name="intersectionPoints"></param>
        /// <param name="pointCount"></param>
        public static void planeBoxIntersection(Plane plane, Vector3 bbMin, Vector3 bbMax, out Vector3[] intersectionPoints)
        {
            int pointCount = 0;
            float t;

            //Vector3[] iPoints = new Vector3[6];
            List<Vector3> iPoints = new List<Vector3>();

            // X - Axis Edges
            Vector3 dir = new Vector3(bbMax.X - bbMin.X, 0.0f, 0.0f);
           
            Vector3 origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMin.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMin.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMin.X, bbMax.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            // Y - Axis Edges
            dir = new Vector3( 0.0f, bbMax.Y - bbMin.Y, 0.0f);

            origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMax.X, bbMin.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMin.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMax.X, bbMin.Y, bbMax.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            // Z - Axis Edges
            dir = new Vector3(0.0f, 0.0f, bbMax.Y - bbMin.Y);

            origin = bbMin;
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMax.X, bbMin.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMin.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            origin = new Vector3(bbMax.X, bbMax.Y, bbMin.Z);
            if (RayPlaneIntersection(origin, dir, plane, out t) && t >= 0.0f && t <= 1.0f)
                iPoints[pointCount++] = origin + dir * t;

            if (pointCount > 0)
                intersectionPoints = SortPointsAndTesselate(iPoints.Take(pointCount).ToArray(), pointCount);
            else
                intersectionPoints = null;

        }

    }
}

*/