Blog Archive

Sunday, June 17, 2012

WPF 3D Plot -- Multivariate Normal

Multivariate Normal distribution is represented by an array of 3d vector associated PDF:

mu=[2 3 5];
sigma=[1,1.5,.5;1.5 3, 0.5;0.5,0.5,2];  %covariance
r=mvnrnd(mu,sigma,100);

p=mvnpdf(r,mu,sigma);
[r,p]  

We can plot mvnrnd in WPF Helix 3D Viewbox using pdf as "Density Material"

<Window>
....
       <ht:HelixViewport3D ItemsSource="{Binding Objects}"  ShowCoordinateSystem="True"  ShowCameraInfo="True"  />

</Window>

 public ObservableCollection<Visual3D> Objects { get; set; }

For some reason I have to build 3D plot by ThreadPool. Calling directly from ctor does not work.
            t = Task.Factory.StartNew(() =>
            {
            }).ContinueWith((t0) =>
            {
                if (Objects.Count == 0)
                    Objects.Add(new DefaultLights());
                Material m = GetData();
                List<Point3DWithDensity> list = Data3D.GetSamplePoints();
                 maxDensity = (from _ in list select _.Density).Max();
                 minDensity = (from _ in list select _.Density).Min();
                foreach( var v in list)
                {
                    BoxVisual3D sp = new BoxVisual3D() { Center = v.Point3D, Height = v.Density, Width = v.Density, Length = v.Density, Material = GetDensityMaterial(v.Density, v.Point3D) };
                 //   SphereVisual3D sp = new SphereVisual3D { Center = v.Point3D, Radius = v.Density, Material = GetDensityMaterial(v.Density, v.Point3D) };
               
                    Objects.Add(sp);
                }
            }, TaskScheduler.FromCurrentSynchronizationContext());


        double maxDensity;
        double minDensity;
        private Material GetDensityMaterial(double d, Point3D p)
        {
            DiffuseMaterial m = new DiffuseMaterial();
            int i = (int)Math.Round(((maxDensity - d) * 1.0 + (d - minDensity) * 10) / (maxDensity - minDensity));
            SolidColorBrush scb = new SolidColorBrush(GetColor(i)); 
            VisualBrush vb = new VisualBrush();
            TextBlock tb = new TextBlock { Text = "(" + p.X.ToString() + "," + p.Y.ToString() + "," + p.Z.ToString() + ")", Background = scb };
           
            vb.Visual = tb;
            m.Brush = vb;
            return m;
        }

    public class Data3D
    {
        public static List<Point3DWithDensity> GetSamplePoints()
        {
            #region Raw Data

            string rawDataFromMATLAB = @"
    2.3631    4.9997    4.8456    0.0129
         ..................
    1.9576    3.5687    1.9473    0.0040
    2.7463    3.6281    4.7240    0.0300";

            #endregion

            List<Point3DWithDensity> list = new List<Point3DWithDensity>();
            string[] pointData = ReplaceMultipleSpace(rawDataFromMATLAB).Split(new string[] { "\r\n" }, StringSplitOptions.RemoveEmptyEntries);
            foreach (string s in pointData)
            {
               string[] xyz= s.Trim().Split(' ');
               list.Add(new Point3DWithDensity()
                {
                    Point3D = new Point3D()
                    {
                        X = Convert.ToDouble(xyz[0].Trim()),
                        Y = Convert.ToDouble(xyz[1].Trim()),
                        Z = Convert.ToDouble(xyz[2].Trim())
                    },
                    Density=Convert.ToDouble(xyz[3].Trim())
                });
            }
            return list;
        }
        static string ReplaceMultipleSpace(string s)
        {
            RegexOptions options = RegexOptions.None; 
            Regex regex = new Regex(@"[ ]{2,}", options);     
            return regex.Replace(s, @" "); 
        }
    }
    public class Point3DWithDensity 
    {
        public Point3D Point3D { get; set; }
        public double Density { get; set; }
    }

        public static Color GetColor(int value) { 
            int startIndex = (value / 10) * 10; 
            int endIndex = startIndex + 10;
            Color startColor = Colors.Red;
            Color endColor = Colors.Blue;
            float weight = (value - startIndex) / (float)(endIndex - startIndex);
          
            return Color.FromArgb(255,
                (byte)Math.Round(startColor.R * (1 - weight) + endColor.R * weight), 
                (byte)Math.Round(startColor.G * (1 - weight) + endColor.G * weight), 
                (byte)Math.Round(startColor.B * (1 - weight) + endColor.B * weight)); 
        } 

No comments: