Data.cs

#region Using Statements
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
using Microsoft.Xna.Framework;
using Microsoft.Xna.Framework.Audio;
using Microsoft.Xna.Framework.Content;
using Microsoft.Xna.Framework.GamerServices;
using Microsoft.Xna.Framework.Graphics;
using Microsoft.Xna.Framework.Input;
using Microsoft.Xna.Framework.Media;
using Microsoft.Xna.Framework.Net;
using Microsoft.Xna.Framework.Storage;
using visLU.Effects;
using Microsoft.Xna.Framework.Graphics.PackedVector;
#endregion

namespace visLU
{
    class Data : VisXnaConrol
    {
        #region variables
        //: change Byte order with NetworkToHostOrder(Int16) if needed

        #region Ray Casting parameters
        public Vector3 stepSize;
        public Vector4 scaleFactor;
        public int iterations;
        public float scaleSteps = 1.0f; //the size of the step to take.  == sample distance
        //Lower means shorter steps. 0.5f or 1.0f are good values.
        public float slice = 0.5f; 
        #endregion

        #region transfer function
        private List<TransferControlPoint> colorControlPoints;
        private List<TransferControlPoint> alphaControlPoints;
        public Texture2D transferTexture;
        public Texture2D transferTextureAlpha; //TDDO:do we use this??? -> if not delete
        #endregion

        #region data
        private List<float> densityValues;
        private float[] densityArray;
        private float[] distinctDensityValuesArray;
        private int[] countDensityValuesArray;

        private Vector3[] gradients;

        public Matrix worldMatrix;
        #endregion


        //declare variables for header information
        ushort sizeX, sizeY, sizeZ;

        private String fileName;

        public Texture3D volumeData;
        public Texture3D volumeGradients; //TDDO: do we use this??? -> if not delete

        public Model cube;

        //private BasicEffect basicEffect; //for debug
        int texturesSaved = 0; //TDDO:do we use this??? -> if not delete

        ContentManager content;

        #region view modes
        private Camera camera;
        private int currentView = 0; //3d view mode
        #endregion

        #endregion

        #region Properties
        public Camera Camera
        {
            set { camera = value; }
            get { return camera; }
        }
        public int CurrentCamera
        {
            set { currentView = value; }
        }

        public Texture3D VolumeData
        {
            get { return volumeData; }
        }
        public Vector3 DataDimensions
        {
            get { return new Vector3(sizeX, sizeY, sizeZ); }
        }
        #endregion

        #region Constructor
        public Data(string _fileName)
        {
            worldMatrix = Matrix.Identity;
            //game.Content.RootDirectory = "Content";
            string a = "" + Directory.GetParent((string)Directory.GetCurrentDirectory());
            string b = "" + Directory.GetParent(a);
            string c = "" + Directory.GetParent(b);

            this.fileName = c + "\\" + _fileName;

            if (string.IsNullOrEmpty(fileName))
            {
                throw new ArgumentNullException(fileName);
            }

            if (!File.Exists(fileName))
            {
                Console.WriteLine(fileName);
                throw new ArgumentException("file not found: " + fileName);
            }

        }
        #endregion

        #region Initialize()
        //-----------------------------------INITIALIZE()-------------------------------------------

        protected override void Initialize()
        {
        }
        //__________________________________________________________________________________________
        #endregion 

        #region LoadContent

        #region LoadContent()
        protected override void LoadContent()
        {
            if (GameProperties.Instance.debugOutput)
            {
                Console.WriteLine("LOADING DATA CONTENT");
            }
            loadDatafromFile();

            #region calculate Ray Casting Parameters
            //calculate stepsize (has to be the ratio of the largest component)
            scaleSteps = GameProperties.Instance.currSampleDist;
            float max = (float)Math.Max(sizeX, Math.Max(sizeY, sizeZ));
            stepSize = Vector3.One / max;

            /*Matrix rotate; -> für debug auskommentriert
            if (GameProperties.Instance.dataFilename.StartsWith("skewed_head"))
            {
                rotate = Matrix.CreateRotationZ(90.0f);
            }*/
            worldMatrix = Matrix.CreateTranslation(new Vector3(-(sizeX / max) / 2.0f,
                                                               -(sizeY / max) / 2.0f,
                                                               -(sizeZ / max) / 2.0f)); //rotate;

            //calculate number of iterations
            iterations = (int)(max * (1.0f / scaleSteps));

            //calculate scaleFactor
            Vector3 sizeVector = new Vector3(sizeX, sizeY, sizeZ);
            Vector3 scale = Vector3.One / ((Vector3.One * max) / sizeVector);
            scaleFactor = new Vector4(scale.X, scale.Y, scale.Z, 1.0f);
            #endregion

            #region content manager
            //get the content manager
            if (GameProperties.Instance.debugOutput)
            {
                Console.WriteLine("creating content manager");
            }
            content = Services.GetService(typeof(ContentManager)) as ContentManager;
            if (content == null)
            {
                //throw new InvalidOperationException("Content Manager Service not found.");
                content = new ContentManager(Services, "Content");
                Services.AddService<ContentManager>(content);
                if (GameProperties.Instance.debugOutput)
                {
                    Console.WriteLine("content manager created");
                }
            }
            try
            {
                cube = content.Load<Model>("Models/box");
                if (GameProperties.Instance.debugOutput)
                {
                    Console.WriteLine("cube loaded");
                }
            }
            catch (Exception e)
            {
                Console.WriteLine("Problem loading cube: " + e);
                throw;
            }
            #endregion

            colorControlPoints = GameProperties.Instance.colorControlPoints;
            alphaControlPoints = GameProperties.Instance.alphaControlPoints;

            transferTexture = new Texture2D(GraphicsDevice, 256, 1, false, SurfaceFormat.Color); ;

            //update slice controls
            setSliceRange(sizeX, sizeY, sizeZ);

        }
        #endregion

        #region loadDatafromFile()
        private void loadDatafromFile()
        {
            if (GameProperties.Instance.debugOutput)
            {
                Console.WriteLine("LOADING DATA FROM FILE");
            }
            #region read data
            Stream file = new FileStream(fileName, FileMode.Open);

            BinaryReader reader = new BinaryReader(file);

            //read header
            sizeX = reader.ReadUInt16();
            sizeY = reader.ReadUInt16();
            sizeZ = reader.ReadUInt16();

            //Console.WriteLine("sizeX =" + sizeX + " sizeY =" + sizeY + " sizeZ =" + sizeZ);

            volumeData = new Texture3D(GraphicsDevice, sizeX, sizeY, sizeZ, false, SurfaceFormat.Rgba64);

            ushort[] buffer = new ushort[sizeX * sizeY * sizeZ];

            int l = 0;
            while (reader.BaseStream.Position < reader.BaseStream.Length)
            {

                buffer[l] = reader.ReadUInt16();
                l++;
            }
            reader.Close();
            #endregion
            if (GameProperties.Instance.debugOutput)
            {
                Console.WriteLine("LOADING DATA FROM FILE");
            }

            #region scale data
            //scale the scalar values to [0, 1]
            densityValues = new List<float>();
            float maxValue = 1.0f / buffer.Max();
            for (int i = 0; i < buffer.Length; i++)
            {
                densityValues.Add((float)buffer[i] * maxValue);
            }
            #endregion

         

            GameProperties.Instance.dataOrigin = new Vector3(0.5f, 0.5f, sizeZ / 2.0f); //TODO: do we need it?

            #region add gradients to VolumeData
            //filter volume data
            if (GameProperties.Instance.enableGaussFilter) { filterGaussian(0.3f); }
            if (GameProperties.Instance.enableAvgFilter) { filterNxNxN(3, true); }


        

            //calculate gradients
            calculateGradients();
            //filter the gradients with an NxNxN box filter
            filterNxNxN(3, false);
            //scale gradients
            scaleGradients();

            #region add black border
            //TODO: better to extend the texture with black border and modify the texture coordinates
            for (int k = 0; k < sizeZ; k++)
            {
                for (int j = 0; j < sizeY; j++)
                {
                    for (int i = 0; i < sizeX; i++)
                    {
                        int idx = i + (j * sizeX) + (k * sizeX * sizeY);
                        if (!((i >= 1 && i < sizeX - 1) &&
                            (j >= 1 && i < sizeY - 1) &&
                            (k >= 1 && k < sizeZ - 1)))
                        {

                            densityValues[idx] = 0.0f;

                        }

                    }
                }
            }
            #endregion

            Rgba64[] densityVectors = new Rgba64[buffer.Length];
            for (int i = 0; i < buffer.Length; i++)
            {
                densityVectors[i] = new Rgba64(gradients[i].X,
                                               gradients[i].Y,
                                               gradients[i].Z,
                                               densityValues[i]);
            }
            #endregion
            //set data to texture

            volumeData.SetData<Rgba64>(densityVectors);

            #region sort data
            densityValues.Sort();
            sortDensityValues();

            /*for (int i = 480000; i < 489000; i++) // for debug
            {
                Console.WriteLine(densityValues[i]);
            }*/
            #endregion
            if (GameProperties.Instance.debugOutput)
            {
                Console.WriteLine("DATA LOADED FROM FILE");
            }
        }
        #endregion

        #region sortDensityValues()
        //: DEBUG THIS
        private void sortDensityValues()
        {

            List<float> distinctDensityValues = new List<float>();
            List<int> countDensityValues = new List<int>();

            #region debug
            /*int countAll = densityValues.Count();
            int testCountAll = 0;
            int totalValues = 0;*/
            #endregion

            for (int i = 0; i < densityValues.Count(); i++)
            {
                int value1 = (int)(densityValues[i] * 255);
                int count = 1;

                while (i < (densityValues.Count() - 1))
                {
                    int value2 = (int)(densityValues[i + 1] * 255);
                    if ((value1 == value2))
                    {
                        count++;
                        i++;
                    }
                    else
                    {
                        break;
                    }

                }

                #region debug
                /*testCountAll += count;
                totalValues++;
                Console.WriteLine("value" + i + " = " + value1 + "   x" + count);*/
                #endregion

                distinctDensityValues.Add(value1);
                countDensityValues.Add(count);

            }
            distinctDensityValuesArray = distinctDensityValues.ToArray();
            countDensityValuesArray = countDensityValues.ToArray();

            GameProperties.Instance.distinctDensityValues = distinctDensityValuesArray;
            GameProperties.Instance.countDensityValues = countDensityValuesArray;

        }
        #endregion

        #endregion

        #region Update()
        protected override void Update() //OLD STUFF->DELETE? -> Leave function empty, it should be overridden by this class 
        {
            #region DEBUG
            //Console.WriteLine("Data Update Call");
            //GameProperties.Instance.countDataUpdates++;
            //Console.WriteLine("Doing Data Update Number " + GameProperties.Instance.countDataUpdates);
            #endregion

            #region move camera TODO
            /*currentMouseState = Mouse.GetState();
            if (currentMouseState.LeftButton == ButtonState.Pressed)
            {
                float fX, fY;
                fX = MathHelper.ToRadians(currentMouseState.X - lastMouseState.X);
                fY = MathHelper.ToRadians(currentMouseState.Y - lastMouseState.Y);

                 *= Matrix.CreateRotationX(fX) * Matrix.CreateRotationY(-fY);

            }
            lastMouseState = currentMouseState;*/
            #endregion

            #region set parameters for rendering
            /*GraphicsDevice.BlendState = BlendState.Additive;
            shader.CurrentTechnique = shader.Techniques["RenderPosition"];
            shader.SetEffectParameter(worldMatrix * camera.ViewProjection, worldMatrix);
            shader.SetEffectParameter(camera.CameraPosition);*/
            #endregion

        }
        #endregion

        #region Draw

        #region Draw()
        protected override void Draw() //OLD STUFF->DELETE? -> Leave function empty, it should be overridden by this class 
        {
            #region DEBUG
            //Console.WriteLine("Draw()");
            //Console.WriteLine("Data Draw Call");
            //Console.WriteLine("3 DRAWING VOLUME for cam=" + currentView + "-----------------------------------------------------");
            #endregion

            DrawCoordinateSystem(); //for debug

        }
        #endregion



        #region DrawCoordinateSystem()
        public void DrawCoordinateSystem()
        {
            Console.WriteLine("DRAWING COORDINATE SYSTEM");
            #region coordinate system
            //Camera camera = (Camera)Game.Services.GetService(typeof(Camera));
            BasicEffect basicEffect = new BasicEffect(GraphicsDevice);
            GraphicsDevice.DepthStencilState = DepthStencilState.None;

            basicEffect.LightingEnabled = false;
            basicEffect.TextureEnabled = false;
            basicEffect.VertexColorEnabled = true;
            basicEffect.World = Matrix.Identity;
            basicEffect.Projection = camera.Projection;
            basicEffect.View = camera.View;

            foreach (EffectPass pass in basicEffect.CurrentTechnique.Passes)
            {
                pass.Apply();
                VertexPositionColor[] pointList = { 
                                      //Z axis
                                      new VertexPositionColor(new Vector3(0.0f, 0.0f, 0.0f), Color.DarkBlue),
                                      new VertexPositionColor(new Vector3(0.0f, 0.0f, -3000.0f), Color.DarkBlue),
                                      new VertexPositionColor(new Vector3(0.0f, 0.0f, 0.0f), Color.Blue), 
                                      new VertexPositionColor(new Vector3(0.0f, 0.0f, 3000.0f), Color.Blue),
                                      //Y axis
                                      new VertexPositionColor(new Vector3(0.0f, 0.0f, 0.0f), Color.DarkGreen),
                                      new VertexPositionColor(new Vector3(0.0f, -3000.0f, 0.0f), Color.DarkGreen),
                                      new VertexPositionColor(new Vector3(0.0f, 0.0f, 0.0f), Color.Green),
                                      new VertexPositionColor(new Vector3(0.0f, 3000.0f, 0.0f), Color.Green),
                                      //X axis
                                      new VertexPositionColor(new Vector3(0.0f, 0.0f, 0.0f), Color.DarkRed),
                                      new VertexPositionColor(new Vector3(-3000.0f, 0.0f, 0.0f), Color.DarkRed),
                                      new VertexPositionColor(new Vector3(0.0f, 0.0f, 0.0f), Color.Red),
                                      new VertexPositionColor(new Vector3(3000.0f, 0.0f, 0.0f), Color.Red),
                                      //X
                                      new VertexPositionColor(new Vector3(10.0f + 2.5f, 5.0f, 0.0f), Color.DarkRed),
                                      new VertexPositionColor(new Vector3(10.0f + 1.25f, 1.25f, 0.0f), Color.DarkRed),
                                      new VertexPositionColor(new Vector3(10.0f + 1.25f, 5.0f, 0.0f), Color.DarkRed),
                                      new VertexPositionColor(new Vector3(10.0f + 2.5f, 1.25f, 0.0f), Color.DarkRed),
                                      //Y
                                      new VertexPositionColor(new Vector3(1.875f, 10.0f + 3.125f, 0.0f), Color.DarkGreen),
                                      new VertexPositionColor(new Vector3(1.875f, 10.0f + 1.125f, 0.0f), Color.DarkGreen),
                                      new VertexPositionColor(new Vector3(1.875f, 10.0f + 3.125f, 0.0f), Color.DarkGreen),
                                      new VertexPositionColor(new Vector3(2.5f, 10.0f + 5.0f, 0.0f), Color.DarkGreen),
                                      new VertexPositionColor(new Vector3(1.875f, 10.0f + 3.125f, 0.0f), Color.DarkGreen),
                                      new VertexPositionColor(new Vector3(1.125f, 10 + 5.0f, 0.0f), Color.DarkGreen),
                                      //Z
                                       new VertexPositionColor(new Vector3(0.0f, 5.0f, 10.0f + 2.5f), Color.DarkBlue),
                                      new VertexPositionColor(new Vector3(0.0f, 5.0f, 10.0f + 1.25f), Color.DarkBlue),
                                      new VertexPositionColor(new Vector3(0.0f, 1.25f, 10.0f + 2.5f), Color.DarkBlue),
                                      new VertexPositionColor(new Vector3(0.0f, 1.25f, 10.0f + 1.25f), Color.DarkBlue),
                                      new VertexPositionColor(new Vector3(0.0f, 1.25f, 10.0f + 2.5f), Color.DarkBlue),
                                      new VertexPositionColor(new Vector3(0.0f, 5.0f, 10.0f + 1.25f), Color.DarkBlue)
                                     
                                      };
                short[] lineListIndices = new short[28] { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27 };
                GraphicsDevice.DrawUserIndexedPrimitives<VertexPositionColor>(
                    PrimitiveType.LineList,
                    pointList,
                    0,  // vertex buffer offset to add to each element of the index buffer
                    28,  // number of vertices in pointList
                    lineListIndices,  // the index buffer
                    0,  // first index element to read
                    14   // number of primitives to draw
                );
            }
            GraphicsDevice.DepthStencilState = DepthStencilState.Default;
            #endregion
        }
        #endregion

        #endregion

        #region transfer function
        public void setControlpoints()
        {
            colorControlPoints = GameProperties.Instance.colorControlPoints;
            alphaControlPoints = GameProperties.Instance.alphaControlPoints;

            float amount = 0.0f;
            Color[] colorArray = new Color[256];
            float[] alphaArray = new float[256];

            #region interpolate between ControlPoints

            // sort ControlPoints by isoValue
            colorControlPoints.Sort(delegate(TransferControlPoint p1,
                TransferControlPoint p2) { return p1.isoValue.CompareTo(p2.isoValue); });
            alphaControlPoints.Sort(delegate(TransferControlPoint p1,
                TransferControlPoint p2) { return p1.isoValue.CompareTo(p2.isoValue); });

            // before first ControlPoint: fill with black and alpha 0.0f
            for (int i = 0; i < colorControlPoints[0].isoValue; i++)
            {
                colorArray[i] = new Color(0.0f, 0.0f, 0.0f, 0.0f);
                alphaArray[i] = 0.0f;
            }

            for (int j = 0; j < colorControlPoints.Count - 1; j++)
            {
                for (int i = colorControlPoints[j].isoValue; i <= colorControlPoints[j + 1].isoValue; i++)
                {
                    amount = (((float)i - colorControlPoints[j].isoValue) /
                        (colorControlPoints[j + 1].isoValue - colorControlPoints[j].isoValue));
                    colorArray[i] = Color.Lerp(new Color(colorControlPoints[j].color),
                    new Color(colorControlPoints[j + 1].color), amount);

                    //interpolate alpha
                    alphaArray[i] = alphaControlPoints[j].color.W * (1.0f - amount) +
                        alphaControlPoints[j + 1].color.W * amount;

                    colorArray[i].A = (byte)(alphaArray[i] * 255);
                }
            }

            //after last ControlPoint: fill with black and alpha
            for (int i = colorControlPoints[colorControlPoints.Count - 1].isoValue + 1; i < 256; i++)
            {
                colorArray[i] = new Color(0.0f, 0.0f, 0.0f, 0.0f);
                alphaArray[i] = 0.0f;
            }
            #endregion

            #region debug
            /*Console.WriteLine("Call Data TF- " + "TF Count:" + colorControlPoints.Count);
            for (int i = 0; i < colorControlPoints.Count; i++)
            {
                Console.WriteLine("CP ISO:" + i + " =" + colorControlPoints[i].isoValue);
                Console.WriteLine("CP COLOR:" + i + " =" + colorControlPoints[i].color);
                Console.WriteLine("CP ALPHA:" + i + " =" + alphaControlPoints[i].color);
                Console.WriteLine("TF for ISO: " + colorControlPoints[i].isoValue + ": " + colorArray[colorControlPoints[i].isoValue]);
            }*/
            #endregion

            #region write to textures
            //write to a textures
            transferTexture = new Texture2D(GraphicsDevice, 256, 1, false, SurfaceFormat.Color);
            transferTextureAlpha = new Texture2D(GraphicsDevice, 256, 1, false, SurfaceFormat.Color);


            Color[] alphaValues = new Color[256];

            for (int i = 0; i < 256; i++)
            {
               
                /*if (i < 5 || i > 250)
                {
                    colorArray[i] = new Color(0.0f, 0.0f, 0.0f, 0.0f);
                    alphaValues[i] = new Color(0.0f, 0.0f, 0.0f, 0.0f);
                    continue;
                }*/
                alphaValues[i] = new Color(alphaArray[i], alphaArray[i], alphaArray[i], alphaArray[i]);

                
            }

            transferTexture.SetData<Color>(colorArray);
            transferTextureAlpha.SetData<Color>(alphaValues);
            //saveScreenshot(transferTexture, "ReferenceTF1", "png");
            #endregion

        }

        #endregion

        #region Gradients

        private float getSampleAt(int x, int y, int z)
        {
            x = (int)MathHelper.Clamp(x, 0, sizeX - 1);
            y = (int)MathHelper.Clamp(y, 0, sizeY - 1);
            z = (int)MathHelper.Clamp(z, 0, sizeZ - 1);

            //y*width -> move rows down
            //z* (width*height) -> move planes into the volume
            int arrayIndex = x + y * sizeX + z * sizeX * sizeY;

            //Console.WriteLine("Sampe at pos:" + arrayIndex + ":" + densityValues[arrayIndex]);
            
            return densityValues[arrayIndex];

            //return (float)mScalars[x + (y * mWidth) + (z * mWidth * mHeight)].ToVector4().W;
        }

        private void calculateGradients()
        {
            /*use Central difference in x-, y- & z-direction (in voxel):
             *                           ( f(x+1) - f(x-1) )
             * gradient f(x,y,z) = 1/2 * | f(y+1) - f(y-1) |
             *                           ( f(z+1) - f(z-1) ) 
             * */

            int sampleSize = GameProperties.Instance.gradientSampleSize;

            gradients = new Vector3[sizeX * sizeY * sizeZ];

            Vector3 normal = Vector3.Zero;
            //Vector3 min = Vector3.One;
            //Vector3 max = Vector3.Zero;
            int i = 0;
            for (int z = 0; z < sizeZ; z++)
            {
                for (int y = 0; y < sizeY; y++)
                {
                    for (int x = 0; x < sizeX; x++)
                    {
                        //todo check this!!!!
                        Vector3 gradient = new Vector3(
                           getSampleAt(x + sampleSize, y, z) - getSampleAt(x - sampleSize, y, z),
                           getSampleAt(x, y + sampleSize, z) - getSampleAt(x, y - sampleSize, z),
                           getSampleAt(x, y, z + sampleSize) - getSampleAt(x, y, z - sampleSize));

                        // does not matter how long gradient is => normalize (instead of * 1/2)
                        gradients[i] = Vector3.Normalize(gradient);
                        //gradients[i].X = gradient.Length(); gradients[i].Y = gradient.Length(); gradients[i].Z = gradient.Length();
                        //or use magnitude az normal if results are not ok
                        
                        //scale and bias to [0,1]
                        //gradients[i].X = (gradients[i].X + 1.0f) / 2.0f; gradients[i].Y = (gradients[i].Y + 1.0f) / 2.0f; gradients[i].Z = (gradients[i].Z + 1.0f) / 2.0f;
                        if (float.IsNaN(gradients[i].X)) gradients[i] = Vector3.Zero;
                        //Console.WriteLine("Gradient " + i + ":" + gradients[i].X + ", " + gradients[i].Y + ", " + gradients[i].Z);
                        //min.X = Math.Min(min.X, gradients[i].X); min.Y = Math.Min(min.Y, gradients[i].Y); min.Z = Math.Min(min.Z, gradients[i].Z);
                        //max.X = Math.Max(max.X, gradients[i].X); max.Y = Math.Max(max.Y, gradients[i].Y); max.Z = Math.Max(max.Z, gradients[i].Z);
                        i++;
                    }
                }
            }
        }


        private void scaleGradients()
        {
            for (int i = 0; i < gradients.Length; i++)
            {
                //scale and bias to [0,1]
                gradients[i].X = (gradients[i].X + 1.0f) / 2.0f; gradients[i].Y = (gradients[i].Y + 1.0f) / 2.0f; gradients[i].Z = (gradients[i].Z + 1.0f) / 2.0f;
                      
            }
        }
        private void filterNxNxN(int n = 3, bool filterData = false)
        {
            List<float> valuesTmpSingle = new List<float>();
            List<Vector3> valuesTmp = new List<Vector3>(); 
            int i = 0;
            for (int z = 0; z < sizeZ; z++)
            {
                for (int y = 0; y < sizeY; y++)
                {
                    for (int x = 0; x < sizeX; x++)
                    {
                        if (filterData)
                        {
                            valuesTmpSingle.Add(sampleNxNxNSingle(x, y, z, n));
                        }
                        else
                        {
                            valuesTmp.Add(sampleNxNxN(x, y, z, n));
                        }
                        i++;
                    }
                }
            }
            if (filterData) { densityValues = valuesTmpSingle; }
            else {gradients = valuesTmp.ToArray<Vector3>(); }
        }

        private Vector3 sampleNxNxN(int x, int y, int z, int n)
        {
            n = (n - 1) / 2;

            Vector3 average = Vector3.Zero;
            int num = 0;

            for (int k = z - n; k <= z + n; k++)
            {
                for (int j = y - n; j <= y + n; j++)
                {
                    for (int i = x - n; i <= x + n; i++)
                    {
                        if (isInBounds(i, j, k))
                        {
                            average += gradients[i + (j * sizeX) + (k * sizeX * sizeY)];
                            num++;
                        }
                    }
                }
            }

            average /= (float)num;
            if (average.X != 0.0f && average.Y != 0.0f && average.Z != 0.0f)
                average.Normalize();

            return average;
        }

        private bool isInBounds(int x, int y, int z)
        {
            return ((x >= 0 && x < sizeX) &&
                    (y >= 0 && y < sizeY) &&
                    (z >= 0 && z < sizeZ));
        }

        private float sampleNxNxNSingle(int x, int y, int z, int n)
        {
            n = (n - 1) / 2;

            float average = 0.0f;
            int num = 0;

            for (int k = z - n; k <= z + n; k++)
            {
                for (int j = y - n; j <= y + n; j++)
                {
                    for (int i = x - n; i <= x + n; i++)
                    {
                        if (isInBounds(i, j, k))
                        {
                            average += densityValues[i + (j * sizeX) + (k * sizeX * sizeY)];
                            num++;
                        }
                    }
                }
            }

            average /= (float)num;
            return average;
        }

        private void filterGaussian(float sigma)
        {
            List<float> valuesTmpSingle = new List<float>();
            float[] gaussianKernel = calculateGaussianKernel(sigma);
            int i = 0;
            for (int z = 0; z < sizeZ; z++)
            {
                for (int y = 0; y < sizeY; y++)
                {
                    for (int x = 0; x < sizeX; x++)
                    {

                        valuesTmpSingle.Add(sampleGaussKernel(x, y, z, gaussianKernel));
                        
                        i++;
                    }
                }
            }
            densityValues = valuesTmpSingle;
        }

        private float sampleGaussKernel(int x, int y, int z, float[] gaussianKernel)
        {
            int n = (int)Math.Pow(gaussianKernel.Length, 1 / 3);
            n = (n - 1) / 2;

            float average = 0.0f;
            int kk = 0; int jk = 0; int ik = 0;

            for (int k = z - n; k <= z + n; k++)
            {
                for (int j = y - n; j <= y + n; j++)
                {
                    for (int i = x - n; i <= x + n; i++)
                    {
                        if (isInBounds(i, j, k))
                        {
                            average += densityValues[i + (j * sizeX) + (k * sizeX * sizeY)] * gaussianKernel[ik + (jk * n) + (kk * n * n)];
                           
                        }
                        ik++;
                    }
                    jk++;
                }
                kk++;
            }

            return average;
        }


        private float[] calculateGaussianKernel(float sigma)
        {
          //obtain minimum kernel size
          int kernelSize = (int)(6 * (sigma) - 1);
          if (kernelSize < 3) kernelSize = 3;

          int width = kernelSize; int height = kernelSize; int depth = kernelSize;

          float[] gaussianKernel = new float[width * height * depth];

          int x_offset = (-1) * (width - 1) / 2;
          int y_offset = (-1) * (height - 1) / 2;
          int z_offset = (-1) * (depth - 1) / 2;

          float norm = 0.0f;

          for (int k = 0; k < depth; k++)
          {
              y_offset = (-1) * (width - 1) / 2;
              for (int j = 0; j < width; j++)
              {
                  x_offset = (-1) * (width - 1) / 2;
                  for (int i = 0; i < width; i++)
                  {
                      int idx = k * (width * width) + j * (width) + i;

                      //obtain kernel coeff
                      gaussianKernel[idx] = 1 / (float)(Math.Sqrt(2 * Math.PI) * sigma * Math.Sqrt(2 * Math.PI) * sigma * Math.Sqrt(2 * Math.PI) * sigma * Math.Exp(-((x_offset * x_offset + y_offset * y_offset + z_offset * z_offset) / (2 * sigma * sigma))));

                      norm = norm + gaussianKernel[idx];

                      x_offset = x_offset + 1;
                  }
                  y_offset = y_offset + 1;
              }
              z_offset = z_offset + 1;
          }

          //normalize to one
          for (int i = 0; i < gaussianKernel.Length; i++)
          {
              gaussianKernel[i] = gaussianKernel[i] / norm;
          }

          return gaussianKernel;
        }

        #endregion

        #region help functions
        private void saveScreenshot(Texture2D texture, String name)
        {
            try
            {
                FileStream stream = File.OpenWrite(name + ".png");
                texture.SaveAsPng(stream, texture.Width, texture.Height);

            }
            catch (Exception e)
            {
                Console.WriteLine("Save Screenshot: ");
                Console.WriteLine(e);
            }
        }

        private void saveScreenshot(Texture2D texture, String name, String ext)
        {
            Console.WriteLine("saving " + name);
            try
            {
                if (!name.ElementAt(name.Length - 4).Equals(".")) name = name + "." + ext;
                FileStream stream = File.OpenWrite(name);
                switch (ext)
                {
                    case "png":
                        texture.SaveAsPng(stream, texture.Width, texture.Height);
                        break;
                    case "jpg":
                        texture.SaveAsJpeg(stream, texture.Width, texture.Height);
                        break;
                    default:
                        if (name.ElementAt(name.Length - 4).Equals(".")) name = name.Substring(0, name.Length - 4);
                        name = name + ".png";
                        texture.SaveAsPng(stream, texture.Width, texture.Height);
                        break;
                }

            }
            catch (Exception e)
            {
                Console.WriteLine("Save Screenshot: ");
                Console.WriteLine(e);
            }

        }
       

        private void setSliceRange(int xRange, int yRange, int zRange)
        {
            GameProperties.Instance.xDataRange = xRange;
            GameProperties.Instance.yDataRange = yRange;
            GameProperties.Instance.zDataRange = zRange;
            //GameProperties.Instance.updateGuiForm = true; old stuff -> delete TODO: do we need this? if not delete
        }
        #endregion
    }
}