Mean Shift

적재된 화물들을 POD(양하지), 화물 타입별로 그룹핑하고, 밀도? 분포에 따라 한덩어리씩 묶어서 그려주는게 목적

역시 C#에는 이런게 없어서 Python을 번역한다.

https://spin.atomicobject.com/2015/05/26/mean-shift-clustering/

https://github.com/mattnedrich/MeanShift_py

아래는 코드전체

namespace MeanShift

{

    public class MeanShiftResult

    {

        public List<PointF> originalPoints;

        public List<PointF> shiftedPoints;

        public List<int> clusterIds;

        public MeanShiftResult(List<PointF> originalPoints, List<PointF> shiftedPoints, List<int> clusterIds)

        {

            this.originalPoints = originalPoints;

            this.shiftedPoints = shiftedPoints;

            this.clusterIds = clusterIds;

        }

    }

    public class MeanShift

    {

        private const double MIN_DISTANCE = 0.0001;

        private const double GROUP_DISTANCE_TOLERANCE = 0.1;

        bool isMultivariateGaussian;

        public MeanShift(bool isMultivariateGaussian = false)

        {

            this.isMultivariateGaussian = isMultivariateGaussian;

        }

        public MeanShiftResult Cluster(List<PointF> points, float kernelBandwidth = 3)

        {

            float maxMinDist = 1;

            int iterationNumber = 0;

            List<PointF> shiftPoints = new List<PointF>();

            foreach(PointF p in points)

                shiftPoints.Add(p);

            List<bool> stillShifting = new List<bool>();

            for (int i = 0; i < points.Count; i++)

                stillShifting.Add(true);

            while(maxMinDist > MIN_DISTANCE){

                maxMinDist = 0;

                iterationNumber += 1;

                for (int i = 0; i < shiftPoints.Count; i++)

                {

                    if(!stillShifting[i])

                        continue;

                    PointF pNew = shiftPoints[i];

                    PointF pNewStart = pNew;

                    pNew = ShiftPoint(pNew, points, kernelBandwidth);

                    float dist = EuclideanDist(pNew, pNewStart);

                    if (dist > maxMinDist)

                        maxMinDist = dist;

                    if(dist < MIN_DISTANCE)

                        stillShifting[i] = false;

                    shiftPoints[i] = pNew;

                }

            }

            List<int> groupAssignments = GroupPoints(shiftPoints);

            return new MeanShiftResult(points, shiftPoints, groupAssignments);

        }

        private PointF ShiftPoint(PointF point,List<PointF> points, float kernelBandwidth)

        {

            float shiftX = 0;

            float shiftY = 0;

            float scaleFactor = 0;

            foreach(PointF pTemp in points)

            {

                float dist = EuclideanDist(point, pTemp);

                float weight = GaussianKernel(dist, kernelBandwidth);

                shiftX += pTemp.X * weight;

                shiftY += pTemp.Y * weight;

                scaleFactor += weight;

            }

            shiftX = shiftX / scaleFactor;

            shiftY = shiftY / scaleFactor;

            return new PointF(shiftX, shiftY);

        }

        

        #region Utils

        private float EuclideanDist(PointF pNew,PointF pNewStart)

        {

            double total = 0;

            total += Math.Pow(pNew.X - pNewStart.X, 2);

            total += Math.Pow(pNew.Y - pNewStart.Y, 2);

            return Convert.ToSingle(Math.Sqrt(total));

        }

        private float GaussianKernel(float distance, float bandwidth)

        {

            double euclideanDistance = Math.Sqrt(Math.Pow(distance, 2)); // 뭔가 빠졌다

            double val = (1/(bandwidth*Math.Sqrt(2*Math.PI))) * Math.Exp(-0.5*Math.Pow(((euclideanDistance/bandwidth)), 2));

            return Convert.ToSingle(val);

        }

        //private float MultivariateGaussianKernel(List<float> distances, List<float> bandwidths)

        //{

        //    //Number of dimensions of the multivariate gaussian

        //    int dim = bandwidths.Count;

        //    //Covariance matrix

        //    double cov = np.multiply(np.power(bandwidths, 2), np.eye(dim))

        //    //Compute Multivariate gaussian (vectorized implementation)

        //    exponent = -0.5 * np.sum(np.multiply(np.dot(distances, np.linalg.inv(cov)), distances), axis=1)

        //    val = (1 / np.power((2 * math.pi), (dim/2)) * np.power(np.linalg.det(cov), 0.5)) * np.exp(exponent)

        //    return val;

        //}

        #endregion

        #region Point Grouper

        private List<int> GroupPoints(List<PointF> points)

        {

            List<int> groupAssignment = new List<int>();

            List<List<PointF>> groups = new List<List<PointF>>();

            int groupIndex = 0;

            foreach(PointF point in points){

                int nearestGroupIndex = DetermineNearestGroup(point, groups);

                if(nearestGroupIndex  == -1)

                {

                    //Create New Group

                    groups.Add(new List<PointF>(){point});

                    groupAssignment.Add(groupIndex);

                    groupIndex += 1;

                }

                else

                {

                    groupAssignment.Add(nearestGroupIndex);

                    groups[nearestGroupIndex].Add(point);

                }

            }

            return groupAssignment;

        }

        private int DetermineNearestGroup(PointF point, List<List<PointF>> groups)

        {

            int nearestGroupIndex = -1;

            int index = 0;

            foreach(List<PointF> group in groups)

            {

                float distanceToGroup = DistanceToGroup(point, group);

                if(distanceToGroup < GROUP_DISTANCE_TOLERANCE)

                    nearestGroupIndex = index;

                index += 1;

            }

            return nearestGroupIndex;

        }

        

        private float DistanceToGroup(PointF point,List<PointF> group)

        {

            float minDistance = float.MaxValue;

            foreach(PointF p in group){

                float dist = EuclideanDist(point, p);

                if(dist < minDistance)

                    minDistance = dist;

            }

            return minDistance;

        }

        #endregion

    }

}