KernelDensityEstimator.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace EdgeClustering
{
    /// \author Manuel Keglevic
    /// \author Thomas Schulz
    /// <summary>
    /// This utility class provides functions for calculating a gaussian kernel and 
    /// a Kernel Density Estimator used to determine the main direction in a specified segment.
    /// </summary>
    static class KernelDensityEstimator
    {
        private static int histogramSize = 360;
        public static int HistogramSize
        {
            get { return KernelDensityEstimator.histogramSize; }
            set { KernelDensityEstimator.histogramSize = value; }
        }

        private static double[] kernel;

        /// <summary>
        /// Generates a Gaussian Kernel with the size depending on 
        /// the given radius and the fixed histogram size
        /// </summary>
        /// <param name="radius">radius in degrees for the gaussian kernel</param>
        static public void GenerateKernel(double radius)
        {
            double h = 180.0 / (double)histogramSize;

            int width = (int)Math.Floor(radius/h) * 2 + 1;
            double sigma = (double)radius / 3.0;
            double norm = 1.0 / (Math.Sqrt(2.0 * Math.PI) * sigma);
            double coeff = 2.0 * sigma * sigma;

            kernel = new double[width];

            double total = 0;
            double x = 0;

            for (int i = 0; i < width; i++)
            {
                x = -radius + i*2*radius/(width-1);
              
                kernel[i] = norm * Math.Exp(-x * x / coeff);
                total += kernel[i];
            }

            for (int i = 0; i < width; i++)
                kernel[i] /= total;
        }
        /// <summary>
        /// Calculates the maximum of the Kernel Density Estimator for the given angles.
        /// </summary>
        /// <param name="angle">List of angles between 0 (inclusive) and 180 (exclusive) degrees.</param>
        public static double CalculateMaximum(double[] angles, double threshold)
        {
            //Calculate histogram
            double[] hist = new double[histogramSize];

            foreach (double a in angles)
            {
                
                //
                //K passt auf keine Fall
                int k = (int)Math.Floor(a * histogramSize/180.0);
                hist[k]++;
            }

            double[] axis = new double[hist.Length];

            int radius = (int)Math.Floor(kernel.Length/2.0);
            double peakValue = 0;
            int peak = 0;

            //Glätten des Histogramms
            for (int j=0; j<hist.Length; j++)   //Loop for each histogram element
            {
                if(hist[j] > 0)
                    for (int i = -radius; i <= radius; i++ )    //Loop for each kernel element
                    {
                        int k = j + i;

                        if (k < 0)
                            k = axis.Length + k;
                        if (k >= axis.Length)
                            k = k - axis.Length;

                        axis[k] += kernel[i + radius] * hist[j];
                       
                        //Calculate Peak
                        if (axis[k] > peakValue)
                        {
                            peakValue = axis[k];
                            peak = k;
                        }
                    }
            }

            //Normalize
            peakValue /= angles.Length;

            //Return peak
            return (peakValue > threshold) ? (peak + 0.5) * 180.0 / histogramSize : -1;
        }
    }
}

---