Normally Distributed Random Numbers

So you have a random generator that produces uniformly distributed random numbers from 0 to 1 and you want to create a normal distributed random number generator.

        /// <summary>
        /// Returns a normally distributed random number
        /// </summary>
        /// <param name="mean">The mean value of the distribution</param>
        /// <param name="stdDev">The standard deviation of the distribution</param>
        /// <returns>A normally distributed random number</returns>
        public Double NextDouble(double mean, double stdDev)
        {
            Random rand = new Random();
            double u1 = rand.NextDouble(), u2 = rand.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
                Math.Sin(2.0 * Math.PI * u2);   // random normal(0,1)
            double randNormal =
                mean + stdDev * randStdNormal;  // random normal(mean,stdDev^2)
            return randNormal;
        }

The only problem is that if you want to create lots of random numbers all the logarithms and square roots just bog down the machine especially if you want real time performance. So what do you do ?

Simple, trade memory for computation time. I created a simple class that is used like the Random class of the .NET framework but it internally caches normally distributed numbers upon initialization.

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

namespace Emgu_Tracker_Advanced_V2.Source
{
    public class RandomNormalCached : Random
    {
        #region Variables
        protected int noCached = 101;
        /// <summary>
        /// Array of precomputed number with 0 mean and 1 variance
        /// </summary>
        protected double[] precomputedNumbers;
        /// <summary>
        /// Internal state counter
        /// </summary>
        protected int counter;
        #endregion

        #region Constructor
        public RandomNormalCached()
            : base()
        {
            counter = 0;
            precomputedNumbers = new double[noCached];
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }

        /// <summary>
        /// Constructor with custom number of cached random numbers
        /// </summary>
        /// <param name="_noCached">Number of cached random numbers</param>
        /// <param name="_seed">Random seed initializer</param>
        public RandomNormalCached(int _noCached, int _seed)
            : base(_seed)
        {
            counter = 0;
            noCached = _noCached;
            precomputedNumbers = new double[noCached];
            // compute a number of normal distributed number with zero mean , 
            // and deviation equal to one
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }

        public RandomNormalCached(int _seed)
            : base(_seed)
        {
            counter = 0;
            precomputedNumbers = new double[noCached];
            // compute a number of normal distributed number with zero mean , 
            // and deviation equal to one
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }
        #endregion

        #region Functions
        public Double NextDouble(double mean, double stdDev)
        {
            double result = mean + stdDev * precomputedNumbers[counter];
            counter = (counter + 1) % noCached;
            return result;
        }

        /// <summary>
        /// Returns a nornal distributed random number
        /// </summary>
        /// <param name="mean">The mean/expected value</param>
        /// <param name="stdDev">The standard deviation</param>
        /// <returns>Normal distributed normal number</returns>
        protected Double ComputeNormalDistribution(double mean, double stdDev)
        {
            double u1 = base.NextDouble(), u2 = base.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
                Math.Sin(2.0 * Math.PI * u2);   // random normal(0,1)
            double randNormal =
                mean + stdDev * randStdNormal;  // random normal(mean,stdDev^2)
            return randNormal;
        }
        #endregion
    }
}

Hand Tracking Testing

I finished the Particle filter , I tweaked at least a dozen parameters, so I just changed the detector and tried to track my hand. It still needs a lot of work but it's going well.

PS, haar cascade training failed after 2 weeks of computations, yippie !!!

Particle Filter Code Snippet

Below I present the very core of my particle filter. I use the standard CONDENSATION algorithm with histogram distance for weighting my samples.

 

        /// <summary>
        /// Particle filter prediction step.
        /// Each particle represents a hypothesis. The particle's state is propagated
        /// according to a linear model , it is then diffused using random noise, the new
        /// state is weighted according to histogram distance.
        /// The new states are finanlly sampled and the process repeats
        /// </summary>
        /// <param name="_imageChannels">The current frames channels to be used to update weight</param>
        /// <param name="_histOriginal">The original normalized histogram</param>
        public void PredictionUpdate(Image<Gray, byte>[] _imageChannels, DenseHistogram _histOriginal)
        {
            Particle particle;
            Matrix<float> tmpState;
            float sumWeight = 0.0f;
            int noParticlesCurrent = 0, k = 0;

            // sample, move , and diffuse particles
            for (int i = 0; i < noParticles; i++)
            {
                particle = particles[i];
                particle.State = particlesSampling[i].State;
                // create normal distributed random state
                randState.SetRandNormal(new MCvScalar(0), new MCvScalar(2.0f * (1.0f - particle.Weight)));
                // propagate state and diffuse it
                particle.State = particle.kalmanFilter.Predict() + randState;
                // add weight to distribution
                sumWeight += particle.UpdateWeight(_imageChannels, _histOriginal);
            }

            Array.Sort(particles);
            Array.Reverse(particles);

            while (noParticlesCurrent < noParticles && sumWeight > 0.0f)
            {
                // get number of samples
                if ((float)rand.NextDouble() <= particles[k].Weight)
                {
                    particlesSampling[noParticlesCurrent].State = particles[k].State;
                    noParticlesCurrent++;
                }
                k = (k + 1) % noParticles;
            }
            // sum the top 10% of the particles
            tmpState = new Matrix<float>(6, 1);
            for (int i = 0; i < particles.Length / 10; i++)
                tmpState = tmpState + particles[i].State;
            tmpState /= particles.Length / 10;
            // convert particles to a detection rectangle
            detectionRectangle = Particle.StateToRectangle(tmpState);
        }

Particle Filter

Well, the particle filter is finished and it is great...it is like a real time genetic algorithm with a very short selection of features and a stupid evaluation process.

Nevertheless it is done and it is working.

Strangely I don't feel any satisfaction after finishing this. It's a bit weird.
Anyway , maybe I'll post some code later on. But (there is always a but) there are so many parameters to tweak that it is driving me crazy...I will now select the features to follow, I'm looking for illumination, rotation and position invariant features to use, unfortunately such features don't exist or are hugely expensive to compute.

Gaussian 2d

I needed to compute a 2-dimensional Gaussian distribution which is very common when using Gabor filters.
First google result Custom 2D Gauss  provided a quick solution but upon first look the implementation didn't take advantage of any of matlab's features (i.e matrix manipulation) or included functions so it is a bit slow.

I ended rigging it up a bit making it very fast.

One problem I didn't address is adding a scale for the pixels but that can be fixed inside applications by scaling the standard deviation.

The code in Matlab:

function mat = Gaussian2d(...
    gsize, center, sigmax, sigmay, theta, offset, factor)
% GAUSSIAN2d - Produces a 2-dimensional gaussian matrix
%
%   Input
%   gsize   :  Size of the output 'gauss', should be a 1x2 vector
%   center  :  The center position of the gaussian [row col]
%   sigmax  :  Std. dev. in the X direction
%   sigmay  :  Std. dev. in the Y direction
%   theta   :  Rotation in degrees
%   offset  :  Minimum value in output
%   factor  :  Related to maximum value of output, should be 
%                    different from zero  
%
%   Output
%   mat     : probability vector for each angle bin 
%-------------------------------------------------------------------------
%%  Validate input arguments
    if ~exist('gsize','var')
        error('Error : gsize argument not set');
    end
    
    if ~exist('center','var')
        center = gsize / 2;
    end
    
    if ~exist('sigmax','var')
        sigmax = 1;
    end
    
    if ~exist('sigmay','var')
        sigmay = 1;
    end
    
    if ~exist('theta','var')
        theta = 0;
    end
     
    if ~exist('offset','var')
        offset = 0;
    end
    
    if ~exist('factor','var')
        factor = 1;
    end
    theta  = (theta/180)*pi;  
%-------------------------------------------------------------------------
%%  Compute
    [X,Y] = meshgrid(1:gsize(1),1:gsize(2));
    Xc  = X - center(1);
    Yc  = Y - center(2);
    xm  = (Xc)*cos(theta) - (Yc)*sin(theta);
    ym  = (Xc)*sin(theta) + (Yc)*cos(theta);
    u   = (xm / sigmax).^2 + (ym / sigmay).^2;
    mat = offset + factor * exp(-u/2);
end

The code in C# in the OpenCv framework :

        /// Create a rotated 2d Gaussian distribution 
        /// </summary>
        /// <param name="_width">Width of the result image in pixels</param>
        /// <param name="_height">Height of the result image in pixels</param>
        /// <param name="_centerX">X-coordinate of the center of the distribution in pixels</param>
        /// <param name="_centerY">Y-coordinate of the center of the distribution in pixels</param>
        /// <param name="_sigmaX">Standard deviation of X-axis in pixels</param>
        /// <param name="_sigmaY">Standard deviation of Y-axis in pixels</param>
        /// <param name="_theta">Rotation of the distribution in degrees</param>
        /// <param name="_offset">Resulting value offset</param>
        /// <param name="_factor">Resulting value multiplier</param>
        /// <returns></returns>
        public static Image<Gray, float> Gaussian2D(
            int   _width    = 10, 
            int   _height   = 10,
            float _centerX  = 5, 
            float _centerY  = 5,
            float _sigmaX   = 1.0f, 
            float _sigmaY   = 1.0f,
            float _theta    = 0.0f, 
            float _offset   = 0.0f, 
            float _factor   = 1.0f)
        {
            float[, ,] data = new float[_height, _width, 1];
            float xc, yc, xm, ym , tmpX, tmpY, u,
                cosTheta = (float)Math.Cos(_theta / 180.0f * Math.PI), 
                sinTheta = (float)Math.Sin(_theta / 180.0f * Math.PI);
            _sigmaX = 1.0f / _sigmaX;
            _sigmaY = 1.0f / _sigmaY;
            for (int i = 0; i < _width; i++)
            {
                xc = i - _centerX;
                for (int j = 0; j < _height; j++)
                {
                    yc = j - _centerY;
                    xm = xc * cosTheta - yc * sinTheta;
                    ym = xc * sinTheta + yc * cosTheta;
                    tmpX = xm * _sigmaX;
                    tmpY = ym * _sigmaY;
                    u = tmpX * tmpX + tmpY * tmpY;
                    data[j, i, 0] = _offset + _factor * (float)Math.Exp(-u / 2);
                }
            }
            return new Image<Gray,float>(data);
        }

Metrics, Distances and all things pretty

Today I post a self contained block of code that defines many metrics and distances often used in optimization, machine learning and computer vision.

I don't usually show my code in public...

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Drawing;
using System.Runtime.InteropServices;
using System.Text;
using System.Linq;

namespace Emgu_Tracker_Advanced_V2.Source
{
    public class MathExtended
    {
        /// <summary>
        /// Defines the different metrics for distance 
        /// </summary>
        public enum DISTANCE_METRIC { 
            BRAY_CURTIS, CANBERRA, CHESSBOARD, CORRELATION, 
            EUCLIDEAN, NORMALIZED_SQUARED_EUCLIDEAN, COSINE_SIMILARITY,
            MANHATTAN, SQUARED_EUCLIDEAN};

        /// <summary>
        /// Measure the distance of 2 n-th dimensional vectors using any of the given metrics
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <param name="_metric">Distance metric</param>
        /// <returns>The distance between the two points</returns>
        public static float Distance(float[] _vector1, float[] _vector2, DISTANCE_METRIC _metric = DISTANCE_METRIC.EUCLIDEAN)
        {
            float result = -1.0f;
            switch (_metric)
            {
                case DISTANCE_METRIC.BRAY_CURTIS:
                    result = DistanceBrayCurtis(_vector1, _vector2);
                    break;
                case DISTANCE_METRIC.CANBERRA:
                    result = DistanceCanberra(_vector1, _vector2);
                    break;
                case DISTANCE_METRIC.CHESSBOARD:
                    result = DistanceChessboard(_vector1, _vector2);
                    break;
                case DISTANCE_METRIC.EUCLIDEAN:
                    result = DistanceEuclidean(_vector1, _vector2);
                    break;
                case DISTANCE_METRIC.NORMALIZED_SQUARED_EUCLIDEAN:
                    result = DistanceNormalizedSquaredEuclidean(_vector1, _vector2);
                    break;
                case DISTANCE_METRIC.COSINE_SIMILARITY:
                    result = DistanceCosineSimilarity(_vector1, _vector2);
                    break;
                case DISTANCE_METRIC.CORRELATION:
                    result = DistanceCorrelation(_vector1, _vector2);
                    break;
                case DISTANCE_METRIC.MANHATTAN:
                    result = DistanceManhattan(_vector1,_vector2);
                    break;
                case DISTANCE_METRIC.SQUARED_EUCLIDEAN:
                    result = DistanceSquaredEuclidean(_vector1, _vector2);
                    break;
                default:
                    result = -1;
                    break; 
            }
            return result;
        }

        /// <summary>
        /// Euclidean distance between 2 2-dimensional points
        /// </summary>
        /// <param name="_point1">Point 1</param>
        /// <param name="_point2">Point 2</param>
        /// <returns>Euclidean distance</returns>
        public static float Distance(Point _point1, Point _point2)
        {
            float xDistance = _point1.X - _point2.X;
            float yDistance = _point1.Y - _point2.Y;
            return (float)System.Math.Sqrt(xDistance * xDistance + yDistance * yDistance);
        }

        /// <summary>
        /// Euclidean distance between 2 2-dimensional points
        /// </summary>
        /// <param name="_point1">Point 1</param>
        /// <param name="_point2">Point 2</param>
        /// <returns>Euclidean distance</returns>
        public static float Distance(PointF _point1, PointF _point2)
        {
            float xDistance = _point1.X - _point2.X;
            float yDistance = _point1.Y - _point2.Y;
            return (float)System.Math.Sqrt(xDistance * xDistance + yDistance * yDistance);
        }

        /// <summary>
        /// Euclidean length of n-th dimensional vector
        /// </summary>
        /// <param name="_vector">Vector</param>
        /// <returns>Length of vector</returns>
        public static float Length(float[] _vector)
        {
            float sum = 0.0f;
            if (_vector == null) 
                return -1.0f;
            for (int i = 0; i < _vector.Length; i++)
                sum += _vector[i] * _vector[i];
            return (float)Math.Sqrt(sum);
        }

        /// <summary>
        /// Run through a list of points and find the closest one to the input vector
        /// The metric type is  euclidean distance
        /// </summary>
        /// <param name="point">The point to search for</param>
        /// <param name="pointList">The distance metric</param>
        /// <returns>The minimum distance point</returns>
        public static PointF NearestNeighbour(PointF point, PointF[] pointList)
        {
            return pointList.OrderBy(fromPoint => Distance(point, fromPoint)).First();
        }

        /// <summary>
        /// Run through an array of float vectors and find the nearest one 
        /// to the input vector, given the distance metric type
        /// </summary>
        /// <param name="_vector">The vector to search for</param>
        /// <param name="_vectorsArray">The array of vectors to find the minimum distance to</param>
        /// <param name="_metric">The distance metric</param>
        /// <returns>The minimum distance vector</returns>
        public static float[] NearestNeighbour(
            float[] _vector, 
            float[][] _vectorsArray,
            DISTANCE_METRIC _metric = DISTANCE_METRIC.EUCLIDEAN)
        {
            return _vectorsArray.OrderBy(fromVector => Distance(_vector, fromVector, _metric)).First();
        }

        /// <summary>
        /// Run through an array of float vectors and find the nearest k vectors 
        /// to the input vector, given the distance metric type
        /// </summary>
        /// <param name="_vector">The vector to search for</param>
        /// <param name="_vectorsArray">The array of vectors to find the minimum distance to</param>
        /// <param name="_metric">The distance metric</param>
        /// <returns>The nearest neighbours in (distance vector) key pairs</returns>
        public static KeyValuePair<float, float[]>[] KNearestNeighbour(
            float[] _vector,
            float[][] _vectorsArray,
            int k = 1,
            DISTANCE_METRIC _metric = DISTANCE_METRIC.EUCLIDEAN)
        {
            SortedList<float, float[]> sortedDistances = new SortedList<float, float[]>(_vectorsArray.Length);
            for (int i = 0; i < _vectorsArray.Length; i++)
                sortedDistances.Add(Distance(_vector, _vectorsArray[i], _metric), _vectorsArray[i]);
            return sortedDistances.Take(k).ToArray<KeyValuePair<float, float[]>>();
        }

        /// <summary>
        /// Run through an array of float vectors and find the nearest k vectors 
        /// to the input vector, given the distance metric type
        /// </summary>
        /// <param name="_vector">The vector to search for</param>
        /// <param name="_vectorsArray">The array of vectors to find the minimum distance to</param>
        /// <param name="_metric">The distance metric</param>
        /// <returns>The nearest neighbours in (distance index) key pairs</returns>
        public static KeyValuePair<float, int>[] KNearestNeighboursIndex(
            float[] _vector,
            float[][] _vectorsArray,
            int k = 1,
            DISTANCE_METRIC _metric = DISTANCE_METRIC.EUCLIDEAN)
        {
            SortedList<float, int> sortedDistances = new SortedList<float, int>(_vectorsArray.Length);
            for (int i = 0; i < _vectorsArray.Length; i++)
                sortedDistances.Add(Distance(_vector, _vectorsArray[i], _metric), i);
            return sortedDistances.Take(k).ToArray<KeyValuePair<float,int>>();
        }

        /// <summary>
        /// Sum of absolute differences of two n-th dimensional vectors
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>Sum of absolute differences else if error -1</returns>
        public static float DistanceManhattan(float[] _vector1, float[] _vector2)
        {
            float sum = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
                sum += Math.Abs(_vector1[i] - _vector2[i]);
            return sum;
        }

        /// <summary>
        /// Sum of squared differences of two n-th dimensional vectors
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>Sum of squared differences else if error -1</returns>
        public static float DistanceSquaredEuclidean(float[] _vector1, float[] _vector2)
        {
            float sum = 0.0f, tmp = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
            {
                tmp = (_vector1[i] - _vector2[i]);
                sum += tmp * tmp;
            }
            return sum;
        }

        /// <summary>
        /// Cosine similarity is a measure of similarity between 
        /// two n-th dimensional vectors vectors by measuring the cosine of the angle between them
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>Cosine similarity else if error -1</returns>
        public static float DistanceCosineSimilarity(float[] _vector1, float[] _vector2)
        {
            float sum = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
                sum += (_vector1[i] * _vector2[i]);
            return sum / (MathExtended.Length(_vector1) * MathExtended.Length(_vector2));
        }

        /// <summary>
        /// Euclidean distance between 2 n-dimensional vectors
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>Euclidean distance</returns>
        public static float DistanceEuclidean(float[] _vector1, float[] _vector2)
        {
            float tmp, sum = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
            {
                tmp = _vector1[i] - _vector2[i];
                sum += tmp * tmp;
            }
            return (float)Math.Sqrt(sum);
        }

        /// <summary>
        /// Normalized squared Euclidean distance between 2 n-th dimensional vectors
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>Normalized squared Euclidean distance else if error -1</returns>
        public static float DistanceNormalizedSquaredEuclidean(float[] _vector1, float[] _vector2)
        {
            float sum1 = 0.0f, sum2 = 0.0f, sumTotal1 = 0.0f, sumTotal2 = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
            {
                sum1 += _vector1[i];
                sum2 += _vector2[i];
            }
            sum1 /= _vector1.Length;
            sum2 /= _vector1.Length;
            for (int i = 0; i < _vector1.Length; i++)
            {
                sumTotal1 += (float)Math.Pow(_vector1[i] - sum1 - _vector2[i] + sum2, 2);
                sumTotal2 += (float)Math.Pow(_vector1[i] - sum1, 2) + (float)Math.Pow(_vector2[i] - sum2, 2);
            }
            return sumTotal1 / (2.0f * sumTotal2);
        }

        /// <summary>
        /// Gives the correlation distance between 2 n-th dimensional vectors
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>Correlation distance else if error -1</returns>
        public static float DistanceCorrelation(float[] _vector1, float[] _vector2)
        {
            float sum1 = 0.0f, sum2 = 0.0f, sumTotal1 = 0.0f, sumTotal2 = 0.0f, sumTotal3 = 0.0f, tmp1 = 0.0f, tmp2 = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
            {
                sum1 += _vector1[i];
                sum2 += _vector2[i];
            }
            sum1 /= _vector1.Length;
            sum2 /= _vector1.Length;
            for (int i = 0; i < _vector1.Length; i++)
            {
                tmp1 = (_vector1[i] - sum1);
                tmp2 = (_vector2[i] - sum2);
                sumTotal1 += tmp1 * tmp2;
                sumTotal2 += tmp1 * tmp1;
                sumTotal3 += tmp2 * tmp2;
            }
            sumTotal2 = (float)Math.Sqrt(sumTotal2);
            sumTotal3 = (float)Math.Sqrt(sumTotal3);
            return 1.0f - sumTotal1 / (sumTotal2 * sumTotal3);
        }

        /// <summary>
        /// Gives the canberra distance between 2 n-th dimensional vectors
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>Canberra distance else if error -1</returns>
        public static float DistanceCanberra(float[] _vector1, float[] _vector2)
        {
            float sum = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
                sum += (float)Math.Abs(_vector1[i] - _vector2[i]) / 
                    (float)(Math.Abs(_vector1[i]) + Math.Abs(_vector2[i]));
            return sum;
        }

        /// <summary>
        /// Gives the BrayCurtis distance between 2 n-th dimensional vectors
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>BrayCurtis distance else if error -1</returns>
        public static float DistanceBrayCurtis(float[] _vector1, float[] _vector2)
        {
            float sum1 = 0.0f, sum2 = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
            {
                sum1 += (float)Math.Abs(_vector1[i] - _vector2[i]);
                sum2 += (float)Math.Abs(_vector1[i] + _vector2[i]);
            }
            return sum1 / sum2;
        }

        /// <summary>
        /// Gives the Chessboard distance between 2 n-th dimensional vectors
        /// </summary>
        /// <param name="_vector1">Vector 1</param>
        /// <param name="_vector2">Vector 2</param>
        /// <returns>Chessboard distance else if error -1</returns>
        public static float DistanceChessboard(float[] _vector1, float[] _vector2)
        {
            float tmp = 0.0f;
            if (_vector1 == null || _vector2 == null)
                return -1.0f;
            if (_vector1.Length != _vector2.Length)
                return -1.0f;
            for (int i = 0; i < _vector1.Length; i++)
                tmp = (float)(Math.Max(Math.Abs(_vector1[i]-_vector2[2]), tmp));
            return tmp;
        }
    }
}

Brand new day

Well, as the saying goes no rest for the wicked, I'm giving particle filter another chance and I'm not getting up until it's done.
The paper I'm following is Object Tracking with an Adaptive Color-Based Particle Filter.

Matching the detected objects with the tracked objects is a delicate process.
Currently my approach is weighted average sum of :

• Rectangle similarity (location, width, height) (0.4f)
• Percentage of rectangles intersection (0.3f)
• Histogram compare with  Correlation metric (0.3f)

 Maybe it's too sensitive but it's a good start

        /// <summary>
        /// Given the detected and tracked objects, match each detected object to at most one tracked object
        /// If no matches are found return -1 to the entry for the detected object
        /// </summary>
        /// <param name="_detections">An array of detected rectangle</param>
        /// <param name="_trackedObjects">An array of tracked objects</param>
        /// <returns>An array with the tracked object's index that is matched to the detected object</returns>
        public static int[] MatchDetections(
            Rectangle[] _detections, 
            TrackedObject[] _trackedObjects,
            Image<Gray, byte>[] _imageChannels)
        {
            float[] weights = new float[] { 0.4f, 0.3f, 0.3f };
            bool[] searched = new bool[_detections.Length];
            int[] matches = new int[_detections.Length];
            float[,] metrics = new float[_detections.Length, _trackedObjects.Length];
            Rectangle rectDetection, rectTracked;
            DenseHistogram histDetection, histTracked;

            // for each detected object compute its 
            // matching metric with all the other tracked objects
            for (int i = 0; i < _detections.Length; i++)
            {
                matches[i] = -1;
                rectDetection = _detections[i];
                rectDetection.Inflate(-rectDetection.Width / 10, -rectDetection.Height / 10);
                histDetection = OpenCVExtended.HistogramAtLocation(_imageChannels, rectDetection, true);
                for (int j = 0; j < _trackedObjects.Length; j++)
                {
                    histTracked = _trackedObjects[j].Histogram;
                    rectTracked = _trackedObjects[j].TrackWindow;

                    metrics[i, j] =
                        OpenCVExtended.RectangleSimilarity(rectDetection, rectTracked) * weights[0] +
                        Math.Max(OpenCVExtended.IntersectionAreaPercentage(rectDetection, rectTracked), 
                            OpenCVExtended.IntersectionAreaPercentage(rectTracked, rectDetection)) * weights[1] +
                        (1.0f - (float)CvInvoke.cvCompareHist(histTracked, histDetection, HIST_METHOD)) * weights[2];
                }
            }
            
            // Given the metric matrix make the best matches 
            // between the detected and tracked objects
            for (int k = 0; k < _detections.Length; k++)
            {
                int topDetection = -1, topTracked = -1;
                float topMetric = float.MinValue;
                for (int i = 0; i < _detections.Length; i++)
                    for (int j = 0; j < _trackedObjects.Length; j++)
                        if (topMetric < metrics[i, j])
                        {
                            topMetric = metrics[i, j];
                            topDetection = i;
                            topTracked = j;
                        }
                // if a value is selected then set to minimum float  
                // the row and collumn of the selected value
                if (topMetric > float.MinValue)
                {
                    for (int j = 0; j < _trackedObjects.Length; j++)
                        metrics[topDetection, j] = float.MinValue;
                    for (int i = 0; i < _detections.Length; i++)
                        metrics[i, topTracked] = float.MinValue;
                }
                // if topmetric exceeds a given threshold then make the assignment
                if (topDetection != -1 && topMetric > MIN_SIMILARITY)
                    matches[topDetection] = topTracked;
            }
            return matches;
        }

PS Haar training is still going...

Multiple object tracking

No code left behind

I decided to dump the particle filter for a while until I found a way to improve it (performance wise).

I worked for a bit on the application's GUI and interface, making any little changes easier to track.

In the coming days I will try use a video file input instead of camera input. This way I expect to streamline my experiments.

A little video showing my current results. I will post another one soon with multiple tracks

PS. I'm still waiting for the haar training cascade to finish. It's been 10 day now.

YES !!!

Well, I gave it another try and used weighted mix of different metrics and I got pretty descent results. I got rid of the sudden flips of objects. Since I track similar objects their histograms tend to be close using any distance metric. Not using the histogram at all felt wrong so I added it but with a lower weight coefficient.

Frustrating day

Matching the detected objects with the tracked objects was a disaster.
I used Bhattacharya's histogram metric and Cosine Similarity for the rectangle.
I will read a couple of papers before attempting again.

Detect VS Track

Given N detected objects and M tracked objects, match each detected object to at most one tracked object.
If no match is found for a detected object create a new tracked object for it and begin tracking it.

Computer vision has no paradigm to follow

It really is like walking in the dark, with some flashes of brilliance from other fields.

How many clusters are there ?

As most problems in computer vision, particle filter estimates an underlying probability distribution. A usual model for the distribution is the mixture of gaussians.
Finding the number of modes is solved usually with heuristic methods and that alone promises sub-optimality.

PS I should upload more pictures from my experiments.

Emgu VS OpenCV in C++

The wrapper works great even though is poorly documented.

Everything in OpenCV is still available and I have the luxury of C#.Net and Windows Forms. Examples found online can be ported easily.

OpenCV code in C++ is readable but upon first look looks like a mess.
While in prototype mode I will continue using the Emgu wrapper. 

Particle Filter

I just finished coding a simple implementation of the CONDENSATION algorithm.
The multi hypothesis framework of the PF is much more robust than the Kalman Filter.
I use Mean shift to track the particles around the HSV feature space.
Hopefully I'll find a clever heuristic way to automatically tweak the parameters.

The waiting is killing me

On October 29th I started training an hand detector using the Viola/Jones method and the openCV tools. I used 12k positive and 20k negative samples.
Today the training completed and the results were a bit off. I had a lot of false negatives. I collected lots of texture samples usually found in 3d modelling sites and started the training again with 32k positive and 32k negative samples. Hopefully the results will be better.

In the meantime I will keep using the old detector and concentrate on developing the tracker.

First Post

I intend to use this blog as a personal calendar for my computer projects.