image thumbnail

Mean square displacement analysis of particles trajectories

by

 

08 Mar 2013 (Updated )

A MATLAB class for the mean square displacement analysis of particle trajectories, with a tutorial.

msdanalyzer.m
classdef msdanalyzer
    %%MSDANALYZER a class for simple mean square displacement analysis.
    %
    % msdanalyzer is a MATLAB per-value class dedicated to the mean square
    % displacement analysis of single particle trajectories. After the
    % instantiation, the user can feed the object with the trajectories of
    % several particles. The class allows him to then derive the MSD
    % curves, the velocity autocorrelation functions, to plot them and fit
    % them, and possible correct for drift if required.
    %
    % If you use this tool for your work, we kindly ask you to cite the
    % following article for which it was created:
    % 
    % Tarantino et al. TNF and IL-1 exhibit distinct ubiquitin requirements
    % for inducing NEMO-IKK supramolecular structures. J Cell Biol (2014)
    % vol. 204 (2) pp. 231-45
    %
    % Jean-Yves Tinevez - Institut Pasteur, 2013 - 2014 
    % <tinevez at pasteur dot fr>
    properties (Constant)
        % Tolerance for binning delays together. Two delays will be binned
        % together if they differ in absolute value by less than
        % 10^-TOLERANCE.
        TOLERANCE = 12;
    end
    
    properties (SetAccess = private)
        % The trajectories stored in a cell array, one T x n_dim per particle
        tracks = {};
        % The dimensionality of the problem
        n_dim
        % Spatial units
        space_units;
        % Time units
        time_units
        % Stores the MSD
        msd
        % Stores the velocty correlation
        vcorr
        % Stores the linear fit of MSD vs t
        lfit
        % Stores the linear fit of the log log plot of MSD vs t
        loglogfit
        % Drift movement
        drift
    end
    
    properties (SetAccess = private, GetAccess = private, Hidden = true)
        % If false, msd needs to be recomputed
        msd_valid = false;
        % If false, vcorr needs to be recomputed
        vcorr_valid = false
        
    end
    
    
    %% Constructor
    methods
        
        function obj = msdanalyzer(n_dim, space_units, time_units)
            %%MSDANALYZER Builds a new MSD analyzer object.
            % obj = MSDANALYZER(dimensionality, space__units, time_units)
            % builds a new empty msdanalyzer object with the specified
            % dimensionality (2 for 2D, 3 for 3D, etc...), spatial units
            % and time units. Units are strings, used only for display and
            % plotting.
            
            if ~isreal(n_dim) || ~isscalar(n_dim) || n_dim < 1 || ~(n_dim == floor(double(n_dim)))
                error('msdanalyzer:BadDimensionality', ...
                    'Dimensionality must be a positive integer, got %f.', ...
                    n_dim);
            end
            obj.n_dim = n_dim;
            obj.space_units = space_units;
            obj.time_units = time_units;
        end
    end
    
    %% Public methods
    methods
        
        function obj = addAll(obj, tracks)
            %%ADDALL Add specified trajectories to msd analyzer.
            %
            % obj = obj.addAll(tracks) adds the given tracks the
            % msdanalyzer object obj.
            %
            % Tracks must be specified as a cell array, one array per
            % track. Each track array must be of size N x (Ndim+1) where N
            % is the number of individual measurements of the particle
            % position and Ndim is the dimensionality specified during
            % object construction. Track array must be arranged with time
            % and space as follow: [ Ti Xi Yi ... ] etc..
            %
            % Adding new tracks to an existing object invalidates its
            % stored MSD values, and will cause them to be recalculated if
            % needed.
            
            tracks = tracks(:);
            n_tracks = numel(tracks);
            
            % Check dimensionality & validity
            for i = 1 : n_tracks
                track = tracks{i};
                
                % Dimensionality
                track_dim = size(track, 2);
                if track_dim ~= obj.n_dim + 1
                    error('msdanalyzer:addAll:BadDimensionality', ...
                        'Tracks must be of size N x (nDim+1) with [ T0 X0 Y0 ... ] etc...');
                end
                
                % Non 0-interval
                t = track(:,1);
                if any(diff(t) == 0)
                    error('msdanalyzer:addAll:BadTimeVector', ...
                        'Track %d has a time vector with indentical time points.', i);
                end
                
            end
            
            % Add to track collection
            obj.tracks = [
                obj.tracks;
                tracks
                ];
            obj.msd_valid = false;
            obj.vcorr_valid = false;
            
        end
        
        
        function varargout = plotTracks(obj, ha, indices, corrected)
            %%PLOTTRACKS Plot the tracks stored in this object.
            %
            % obj.plotTracks plots the particle trajectories stored in the
            % msdanalyzer object obj in the current axes. This method only
            % works for 2D or 3D problems.
            %
            % obj.plotTracks(ha) plots the trajectories in the axes with
            % handle ha.
            %
            % obj.plotTracks(ha, indices) where indices is a vector, allows
            % to specify the track to be plotted, using their indices.
            % Leave the vector empty to plot all trajectories.
            %
            % obj.plotTracks(ha, indices, corrected) where corrected is a a
            % boolean flag, allows to specify whether the plot should
            % display the trajectories corrected for drift (true) or
            % uncorrected (false). A proper drift vector must be computed
            % prior to setting this flag to true. See
            % msdanalyzer.computeDrift.
            %
            % hps = obj.plotTracks(...) returns the handles to the
            % individual line objects created.
            %
            % [hps, ha] = obj.plotTracks(...) returns also the handle to
            % the axes handle the trajectories are plot in.
            
            if obj.n_dim < 2 || obj.n_dim > 3
                error('msdanalyzer:plotTracks:UnsupportedDimensionality', ...
                    'Can only plot tracks for 2D or 3D problems, got %dD.', obj.n_dim);
            end
            
            if nargin < 2
                ha = gca;
            end
            if nargin < 3 || isempty(indices)
                indices = 1 : numel(obj.tracks);
            end
            if nargin < 4
                corrected = false;
            end
            
            n_tracks = numel(indices);
            colors = jet(n_tracks);
            
            hold(ha, 'on');
            hps = NaN(n_tracks, 1);
            
            if obj.n_dim == 2
                % 2D case
                for i = 1 : n_tracks
                    
                    index = indices(i);
                    track = obj.tracks{index};
                    
                    x = track(:,2);
                    y = track(:,3);
                    
                    if corrected && ~isempty(obj.drift)
                        tdrift = obj.drift(:,1);
                        xdrift = obj.drift(:, 2);
                        ydrift = obj.drift(:, 3);
                        t = track(:,1);
                        [~, index_in_drift_time, ~] = intersect(tdrift, t);
                        % Subtract drift position to track position
                        x = x - xdrift(index_in_drift_time);
                        y = y - ydrift(index_in_drift_time);
                    end
                    
                    hps(i) =  plot(ha, x, y, ...
                        'Color', colors(i,:));
                    
                end
                
            else
                % 3D case
                for i = 1 : n_tracks
                    
                    index = indices(i);
                    track = obj.tracks{index};
                    
                    x = track(:,2);
                    y = track(:,3);
                    z = track(:,4);
                    
                    if corrected && ~isempty(obj.drift)
                        tdrift = obj.drift(:,1);
                        xdrift = obj.drift(:, 2);
                        ydrift = obj.drift(:, 3);
                        zdrift = obj.drift(:, 4);
                        t = track(:,1);
                        [~, index_in_drift_time, ~] = intersect(tdrift, t);
                        % Subtract drift position to track position
                        x = x - xdrift(index_in_drift_time);
                        y = y - ydrift(index_in_drift_time);
                        z = z - zdrift(index_in_drift_time);
                    end
                    
                    hps(i) =  plot3(ha, x, y, z, ...
                        'Color', colors(i,:));
                    
                end
                
            end
            
            % Output
            if nargout > 0
                varargout{1} = hps;
                if nargout > 1
                    varargout{2} = ha;
                end
            end
            
        end
        
        
        function varargout = plotDrift(obj, ha)
            %%PLOTDRIFT Plot drift stored in this object.
            %
            % obj.plotDrift plots the calculated drift position in the
            % current axes.
            %
            % obj.plotDrift(ha) plots the drift position in the axes
            % specified by the handle ha.
            %
            % hp = obj.plotDrift(...) returns the plot handle for the
            % created line.
            %
            % [hp, ha] = obj.plotDrift(...) also returns the axes handle.
            
            if obj.n_dim < 2 || obj.n_dim > 3
                error('msdanalyzer:plotDrift:UnsupportedDimensionality', ...
                    'Can only plot drift for 2D or 3D problems, got %dD.', obj.n_dim);
            end
            
            if isempty(obj.drift)
                if nargout > 0
                    varargout{1} = [];
                    if nargout > 1
                        varargout{2} = [];
                    end
                end
                return
            end
            
            if nargin < 2
                ha = gca;
            end
            
            if obj.n_dim == 2
                hp = plot(ha, obj.drift(:,2), obj.drift(:,3), 'k-', ...
                    'LineWidth', 2);
            else
                hp = plot3(ha, obj.drift(:,2), obj.drift(:,3), obj.drift(:,4), 'k-', ...
                    'LineWidth', 2);
            end
            
            if nargout > 0
                varargout{1} = hp;
                if nargout > 1
                    varargout{2} = ha;
                end
            end
            
        end
        
        
        function varargout = labelPlotTracks(obj, ha)
            %%LABELPLOTTRACKS A convenience method to set the axes labels.
            %
            % obj.labelPlotTracks(ha) sets the axis label of the axes with
            % the specified handle ha. It is meant for axes containing the
            % plot of the particles trajectories and their drift.
            %
            % hl = obj.plotTracks(...) returns the handle to the generated
            % labels.
            
            if obj.n_dim < 2 || obj.n_dim > 3
                error('msdanalyzer:labelPlotTracks:UnsupportedDimensionality', ...
                    'Can only label axis for 2D or 3D problems, got %dD.', obj.n_dim);
            end
            
            if nargin < 2
                ha = gca;
            end
            
            hl = NaN(obj.n_dim, 1);
            hl(1) = xlabel(ha, ['X (' obj.space_units ')'] );
            hl(2) = ylabel(ha, ['Y (' obj.space_units ')'] );
            if obj.n_dim ==3
                hl(3) = zlabel(ha, ['Z (' obj.space_units ')'] );
            end
            axis equal
            if nargout > 0
                varargout{1} = hl;
            end
        end
        
        
        function varargout = labelPlotMSD(obj, ha)
            %%LABELPLOTMSD A convenience method to set the axes labels.
            %
            % obj.labelPlotMSD(ha) sets the axis label of the axes with
            % the specified handle ha. It is meant for axes containing the
            % plot of the mean-square-displacement.
            %
            % hl = obj.plotMSD(...) returns the handle to the generated
            % labels.
            
            if nargin < 2
                ha = gca;
            end
            
            hl = NaN(2, 1);
            hl(1) = xlabel(ha, ['Delay (' obj.time_units ')']);
            hl(2) = ylabel(ha, ['MSD (' obj.space_units '^2)']);
            
            xl = xlim(ha);
            xlim(ha, [0 xl(2)]);
            yl = ylim(ha);
            ylim(ha, [0 yl(2)]);
            box(ha, 'off')
            
            if nargout > 0
                varargout{1} = hl;
            end
        end
        
        
        function varargout = labelPlotVCorr(obj, ha)
            %%LABELPLOTVCORR A convenience method to set the axes labels.
            %
            % obj.labelPlotVCorr(ha) sets the axis label of the axes with
            % the specified handle ha. It is meant for axes containing the
            % plot of the velocity autocorrelation.
            %
            % hl = obj.labelPlotVCorr(...) returns the handle to the generated
            % labels.
            %
            % [hl, hline] = obj.labelPlotVCorr(...) also returns the handle
            % to generate y=0 dashed line.
            
            if nargin < 2
                ha = gca;
            end
            
            hl = NaN(2, 1);
            hl(1) = xlabel(ha, ['Delay (' obj.time_units ')']);
            hl(2) = ylabel(ha, 'Normalized velocity autocorrelation');
            
            
            xl = xlim(ha);
            xlim(ha, [0 xl(2)]);
            box(ha, 'off')
            
            hline = line([0 xl(2)], [0 0], ...
                'Color', 'k', ...
                'LineStyle', '--');
            uistack(hline, 'bottom')
            
            if nargout > 0
                varargout{1} = hl;
                if nargout > 1
                    varargout{2} = hline;
                end
            end
        end
        
        
        function obj = computeMSD(obj, indices)
            %%COMPUTEMSD Compute the mean-squared-displacement for this object.
            %
            % obj = obj.computeMSD computes the MSD for all the tracks stored
            % in this object. If a drift correction was computed prior to this
            % method call, it is used to correct positions before MSD
            % calculation.
            %
            % Results are stored in the msd field of this object as a cell
            % array, one cell per particle. The array is a double array of size
            % N x 4, and is arranged as follow: [dt mean std N ; ...] where dt
            % is the delay for the MSD, mean is the mean MSD value for this
            % delay, std the standard deviation and N the number of points in
            % the average.
            %
            % obj = obj.computeMSD(indices) computes the MSD only for the
            % particles with the specified indices. Use an empty array to take
            % all particles.
            
            if nargin < 2 || isempty(indices)
                indices = 1 : numel(obj.tracks);
            end
            
            n_tracks = numel(indices);
            fprintf('Computing MSD of %d tracks... ', n_tracks);
            
            % First, find all possible delays in time vectors.
            % Time can be arbitrary spaced, with frames missings,
            % non-uniform sampling, etc... so we have to do this clean.
            % We use a certain tolerance to bin delays together
            delays = obj.getAllDelays;
            n_delays = numel(delays);
            
            obj.msd = cell(n_tracks, 1);
            if ~isempty(obj.drift)
                tdrift = obj.drift(:,1);
                xdrift = obj.drift(:, 2:end);
            end
            
            fprintf('%4d/%4d', 0, n_tracks);
            
            for i = 1 : n_tracks
                
                fprintf('\b\b\b\b\b\b\b\b\b%4d/%4d', i, n_tracks);
                
                mean_msd    = zeros(n_delays, 1);
                M2_msd2     = zeros(n_delays, 1);
                n_msd       = zeros(n_delays, 1);
                
                index = indices(i);
                track = obj.tracks{index};
                t = track(:,1);
                t = msdanalyzer.roundn(t, msdanalyzer.TOLERANCE);
                X = track(:, 2:end);
                
                % Determine drift correction
                if ~isempty(obj.drift)
                    % Determine target delay index in bulk
                    [~, index_in_drift_time, index_in_track_time] = intersect(tdrift, t);
                    % Keep only track times that can be corrected.
                    X = X(index_in_track_time, :);
                    t = t(index_in_track_time);
                    % Subtract drift position to track position
                    X = X - xdrift(index_in_drift_time, :);
                    
                end
                
                
                n_detections = size(X, 1);
                
                for j = 1 : n_detections - 1
                    
                    % Delay in physical units
                    dt = t(j+1:end) - t(j);
                    dt = msdanalyzer.roundn(dt, msdanalyzer.TOLERANCE);
                    
                    % Determine target delay index in bulk
                    [~, index_in_all_delays, ~] = intersect(delays, dt);
                    
                    % Square displacement in bulk
                    dX = X(j+1:end,:) - repmat(X(j,:), [(n_detections-j) 1] );
                    dr2 = sum( dX .* dX, 2);
                    
                    % Store for mean computation / Knuth
                    n_msd(index_in_all_delays)     = n_msd(index_in_all_delays) + 1;
                    delta = dr2 - mean_msd(index_in_all_delays);
                    mean_msd(index_in_all_delays) = mean_msd(index_in_all_delays) + delta ./ n_msd(index_in_all_delays);
                    M2_msd2(index_in_all_delays)  = M2_msd2(index_in_all_delays) + delta .* (dr2 - mean_msd(index_in_all_delays));
                end
                
                n_msd(1) = n_detections;
                std_msd = sqrt( M2_msd2 ./ n_msd ) ;
                
                % We replace points for which N=0 by Nan, to later treat
                % then as missing data. Indeed, for each msd cell, all the
                % delays are present. But some tracks might not have all
                % delays
                delay_not_present = n_msd == 0;
                mean_msd( delay_not_present ) = NaN;
                
                obj.msd{index} = [ delays mean_msd std_msd n_msd ];
                
            end
            fprintf('\b\b\b\b\b\b\b\b\bDone.\n')
            
            obj.msd_valid = true;
            
        end
        
        
        function varargout = plotMSD(obj, ha, indices, errorbar)
            %% PLOTMSD Plot the mean square displacement curves.
            %
            % obj.plotMSD plots the MSD curves in the current axes.
            %
            % obj.plotMSD(ha) plots the MSD curves in the axes
            % specified by the handle ha.
            %
            % obj.plotMSD(ha, indices) plots the MSD curves for the
            % particles with the specified indices only. Leave empty to
            % plot MSD for all particles.
            %
            % obj.plotMSD(ha, indices, errorbar), where errorbar is a
            % boolean flag, allows to specify whether the curves should be
            % plotted with error bars (equal to standard deviation). It is
            % false by default.
            %
            % hps =  obj.plotMSD(...) returns the handle array for the
            % lines generated.
            %
            % [hps, ha] =  obj.plotMSD(...) also return the axes handle in
            % which the lines were plotted.
            
            if ~obj.msd_valid
                obj = obj.computeMSD;
            end
            
            if nargin < 2
                ha = gca;
            end
            if nargin < 3 || isempty(indices)
                indices = 1 : numel(obj.msd);
            end
            if nargin < 4
                errorbar = false;
            end
            
            n_spots = numel(indices);
            colors = jet(n_spots);
            
            hold(ha, 'on');
            if errorbar
            else
                hps = NaN(n_spots, 1);
            end
            
            for i = 1 : n_spots
                
                index = indices(i);
                
                msd_spot = obj.msd{index};
                t = msd_spot(:,1);
                m = msd_spot(:,2);
                if errorbar
                    s = msd_spot(:,3);
                    hps(i) = msdanalyzer.errorShade(ha, t, m, s, colors(i,:), true);
                else
                    hps(i) = plot(ha, t, m, 'Color', colors(i,:));
                end
                
            end
            
            obj.labelPlotMSD(ha);
            
            if nargout > 0
                varargout{1} = hps;
                if nargout > 1
                    varargout{2} = ha;
                end
            end
            
            
        end
        
        
        function obj = fitLogLogMSD(obj, clip_factor)
            %%FITLOGLOGMSD Fit the log-log MSD to determine behavior.
            %
            % obj = obj.fitLogLogMSD fits each MSD curve stored in this object
            % in a log-log fashion. If x = log(delays) and y = log(msd) where
            % 'delays' are the delays at which the msd is calculated, then this
            % method fits y = f(x) by a straight line y = alpha * x + gamma, so
            % that we approximate the MSD curves by MSD = gamma * delay^alpha.
            % By default, only the first 25% of each MSD curve is considered
            % for the fit,
            %
            % Results are stored in the 'loglogfit' field of the returned
            % object. It is a structure with 3 fields:
            % - alpha: all the values for the slope of the log-log fit.
            % - gamma: all the values for the value at origin of the log-log fit.
            % - r2fit: the adjusted R2 value as a indicator of the goodness of
            % the fit.
            %
            % obj = obj.fitLogLogMSD(clip_factor) does the fit, taking into
            % account only the first potion of each MSD curve specified by
            % 'clip_factor' (a double between 0 and 1). If the value
            % exceeds 1, then the clip factor is understood to be the
            % maximal number of point to take into account in the fit. By
            % default, it is set to 0.25.
            
            if nargin < 2
                clip_factor = 0.25;
            end
            
            if ~obj.msd_valid
                obj = obj.computeMSD;
            end
            n_spots = numel(obj.msd);
            
            if clip_factor < 1
                fprintf('Fitting %d curves of log(MSD) = f(log(t)), taking only the first %d%% of each curve... ',...
                    n_spots, ceil(100 * clip_factor) )
            else
                fprintf('Fitting %d curves of log(MSD) = f(log(t)), taking only the first %d points of each curve... ',...
                    n_spots, round(clip_factor) )
            end
            
            alpha = NaN(n_spots, 1);
            gamma = NaN(n_spots, 1);
            r2fit = NaN(n_spots, 1);
            ft = fittype('poly1');
            
            fprintf('%4d/%4d', 0, n_spots);
            for i_spot = 1 : n_spots
                
                fprintf('\b\b\b\b\b\b\b\b\b%4d/%4d', i_spot, n_spots);
                
                msd_spot = obj.msd{i_spot};
                
                t = msd_spot(:,1);
                y = msd_spot(:,2);
                w = msd_spot(:,4);
                
                % Clip data
                if clip_factor < 1
                    t_limit = 2 : round(numel(t) * clip_factor);
                else
                    t_limit = 2 : min(1+round(clip_factor), numel(t));
                end
                t = t(t_limit);
                y = y(t_limit);
                w = w(t_limit);
                
                % Thrash bad data
                nonnan = ~isnan(y);
                
                t = t(nonnan);
                y = y(nonnan);
                w = w(nonnan);
                
                if numel(y) < 2
                    continue
                end
                
                xl = log(t);
                yl = log(y);
                
                bad_log =  isinf(xl) | isinf(yl);
                xl(bad_log) = [];
                yl(bad_log) = [];
                w(bad_log) = [];
                
                if numel(xl) < 2
                    continue
                end
                
                [fo, gof] = fit(xl, yl, ft, 'Weights', w);
                
                alpha(i_spot) = fo.p1;
                gamma(i_spot) = exp(fo.p2);
                r2fit(i_spot) = gof.adjrsquare;
                
            end
            fprintf('\b\b\b\b\b\b\b\b\bDone.\n')
            
            obj.loglogfit = struct(...
                'alpha', alpha, ...
                'gamma', gamma, ...
                'r2fit', r2fit);
            
        end
        
        
        function obj = fitMSD(obj, clip_factor)
            %%FITMSD Fit all MSD curves by a linear function.
            %
            % obj = obj.fitMSD fits all MSD curves by a straight line 
            %                      y = a * x + b.
            % The fit is therefore rigorously valid only for purely
            % diffusive behavior.
            %
            % Results are stored in the 'fit' field of the returned
            % object. It is a structure with 2 fields:
            % - a: all the values of the slope of the linear fit.
            % - b: all the values for the intersect of the linear fit.
            % - r2fit: the adjusted R2 value as a indicator of the goodness 
            % of the fit.
            %
            % obj = obj.fitMSD(clip_factor) does the fit, taking into
            % account only the first potion of the average MSD curve
            % specified by 'clip_factor' (a double between 0 and 1). If the
            % value exceeds 1, then the clip factor is understood to be the
            % maximal number of point to take into account in the fit. By
            % default, it is set to 0.25.
            
            
            if nargin < 2
                clip_factor = 0.25;
            end
            
            if ~obj.msd_valid
                obj = obj.computeMSD;
            end
            n_spots = numel(obj.msd);
            
            if clip_factor < 1
                fprintf('Fitting %d curves of MSD = f(t), taking only the first %d%% of each curve... ',...
                    n_spots, ceil(100 * clip_factor) )
            else
                fprintf('Fitting %d curves of MSD = f(t), taking only the first %d points of each curve... ',...
                    n_spots, round(clip_factor) )
            end
            
            a = NaN(n_spots, 1);
            b = NaN(n_spots, 1);
            r2fit = NaN(n_spots, 1);
            ft = fittype('poly1');
            
            fprintf('%4d/%4d', 0, n_spots);
            for i_spot = 1 : n_spots
                
                fprintf('\b\b\b\b\b\b\b\b\b%4d/%4d', i_spot, n_spots);
                
                msd_spot = obj.msd{i_spot};
                
                t = msd_spot(:,1);
                y = msd_spot(:,2);
                w = msd_spot(:,4);
                
                % Clip data, never take the first one dt = 0
                if clip_factor < 1
                    t_limit = 2 : round(numel(t) * clip_factor);
                else
                    t_limit = 2 : min(1+round(clip_factor), numel(t));
                end
                t = t(t_limit);
                y = y(t_limit);
                w = w(t_limit);
                
                % Thrash bad data
                nonnan = ~isnan(y);
                x = t(nonnan);
                y = y(nonnan);
                w = w(nonnan);
                
                if numel(y) < 2
                    continue
                end
                
                [fo, gof] = fit(x, y, ft, 'Weights', w);
                
                a(i_spot) = fo.p1;
                b(i_spot) = fo.p2;
                r2fit(i_spot) = gof.adjrsquare;
                
            end
            fprintf('\b\b\b\b\b\b\b\b\bDone.\n')
            
            obj.lfit = struct(...
                'a', a, ...
                'b', b, ...
                'r2fit', r2fit);
            
        end
        
        
        function  varargout = fitMeanMSD(obj, clip_factor)
            %%FITMEANMSD Fit the weighted averaged MSD by a linear function.
            %
            % obj.fitMeanMSD computes and fits the weighted mean MSD by a
            % straight line y = a * x. The fit is therefore valid only for
            % purely diffusive behavior. Fit results are displayed in the
            % command window.
            %
            % obj.fitMeanMSD(clip_factor) does the fit, taking into account
            % only the first potion of the average MSD curve specified by
            % 'clip_factor' (a double between 0 and 1). If the value
            % exceeds 1, then the clip factor is understood to be the
            % maximal number of point to take into account in the fit. By
            % default, it is set to 0.25.
            %
            % [fo, gof] = obj.fitMeanMSD(...) returns the fit object and the
            % goodness of fit.
            
            if nargin < 2
                clip_factor = 0.25;
            end
            
            if ~obj.msd_valid
                obj = obj.computeMSD;
            end
            
            ft = fittype('poly1');
            mmsd = obj.getMeanMSD;
            
            t = mmsd(:,1);
            y = mmsd(:,2);
            w = 1./mmsd(:,3);
            
            % Clip data, never take the first one dt = 0
            if clip_factor < 1
                t_limit = 2 : round(numel(t) * clip_factor);
            else
                t_limit = 2 : min(1+round(clip_factor), numel(t));
            end
            t = t(t_limit);
            y = y(t_limit);
            w = w(t_limit);
            
            [fo, gof] = fit(t, y, ft, 'Weights', w);
            
            ci = confint(fo);
            str = sprintf([
                'Estimating D through linear weighted fit of the mean MSD curve.\n', ...
                'D = %.3e with 95%% confidence interval [ %.3e - %.3e ].\n', ...
                'Goodness of fit: R = %.3f.' ], ...
                fo.p1/2/obj.n_dim, ci(1)/2/obj.n_dim, ci(2)/2/obj.n_dim, gof.adjrsquare);
            disp(str)
            
            if nargout > 0
                varargout{1} = fo;
                if nargout > 1
                    varargout{2} = gof;
                end
            end
            
        end
        
        
        function msmsd = getMeanMSD(obj, indices)
            %%GETMEANMSD Compute the weighted mean of all MSD curves.
            %
            % msd = obj.getMeanMSD computes and return the weighted mean of all
            % MSD curves stored in this object. All possible delays are first
            % derived, and for each delay, a weighted mean is computed from all
            % the MSD curves stored in this object. Weights are set to be the
            % number of points averaged to generate the mean square
            % displacement value at the given delay. Thus, we give more weight
            % to MSD curves with greater certainty (larger number of elements
            % averaged).
            %
            % Results are returned as a N x 4 double array, and ordered as
            % following: [ dT M STD N ] with:
            % - dT the delay vector
            % - M the weighted mean of MSD for each delay
            % - STD the weighted standard deviation
            % - N the number of degrees of freedom in the weighted mean
            % (see http://en.wikipedia.org/wiki/Weighted_mean)
            %
            % msd = obj.getMeanMSD(indices) only takes into account the MSD
            % curves with the specified indices.
            
            if ~obj.msd_valid
                obj = obj.computeMSD(indices);
            end
            
            if nargin < 2 || isempty(indices)
                indices = 1 : numel(obj.msd);
            end
            
            n_tracks = numel(indices);
            
            % First, collect all possible delays
            all_delays = cell(n_tracks, 1);
            for i = 1 : n_tracks
                index = indices(i);
                all_delays{i} = obj.msd{index}(:,1);
            end
            delays = unique( vertcat( all_delays{:} ) );
            n_delays = numel(delays);
            
            % Collect
            sum_weight          = zeros(n_delays, 1);
            sum_weighted_mean   = zeros(n_delays, 1);
            
            % 1st pass
            for i = 1 : n_tracks
                
                index = indices(i);
                
                t = obj.msd{index}(:,1);
                m = obj.msd{index}(:,2);
                n = obj.msd{index}(:,4);
                
                % Do not tak NaNs
                valid = ~isnan(m);
                t = t(valid);
                m = m(valid);
                n = n(valid);
                
                % Find common indices
                [~, index_in_all_delays, ~] = intersect(delays, t);
                
                % Accumulate
                sum_weight(index_in_all_delays)           = sum_weight(index_in_all_delays)         + n;
                sum_weighted_mean(index_in_all_delays)    = sum_weighted_mean(index_in_all_delays)  + m .* n;
            end
            
            % Compute weighted mean
            mmean = sum_weighted_mean ./ sum_weight;
            
            % 2nd pass: unbiased variance estimator
            sum_weighted_variance = zeros(n_delays, 1);
            sum_square_weight     = zeros(n_delays, 1);
            
            for i = 1 : n_tracks
                
                index = indices(i);
                
                t = obj.msd{index}(:,1);
                m = obj.msd{index}(:,2);
                n = obj.msd{index}(:,4);
                
                % Do not tak NaNs
                valid = ~isnan(m);
                t = t(valid);
                m = m(valid);
                n = n(valid);
                
                % Find common indices
                [~, index_in_all_delays, ~] = intersect(delays, t);
                
                % Accumulate
                sum_weighted_variance(index_in_all_delays)    = sum_weighted_variance(index_in_all_delays)  + n .* (m - mmean(index_in_all_delays)).^2 ;
                sum_square_weight(index_in_all_delays)        = sum_square_weight(index_in_all_delays)      + n.^2;
            end
            
            % Standard deviation
            mstd = sqrt( sum_weight ./ (sum_weight.^2 - sum_square_weight) .* sum_weighted_variance );
            
            % Output [ T mean std Nfreedom ]
            msmsd = [ delays mmean mstd (sum_weight.^2 ./ sum_square_weight) ];
            
        end
        
        
        function varargout = plotMeanMSD(obj, ha, errorbar, indices)
            %%PLOTMEANMSD Plot the weighted mean of the MSD curves.
            %
            % obj,plotMeanMSD computes and plots the weighted of all MSD
            % curves. See msdanalyzer.getMeanMSD.
            %
            % obj,plotMeanMSD(ha) plots the curve in the axes with the
            % specified handle.
            %
            % obj,plotMeanMSD(ha, errorbar) where 'errorbar' is a boolean allow
            % to specify whether to plot the curve with error bars indicating
            % the weighted standard deviation. Default is false.
            %
            % obj,plotMeanMSD(ha, errorbar, indices) computes and plots the
            % mean only fothe MSD curves whose indices are given on the
            % 'indices' array.
            %
            % h = obj,plotMeanMSD(...) returns the handle to the line plotted.
            %
            % [h, ha] = obj,plotMeanMSD(...) also returns the handle of the
            % axes in which the curve was plotted.
            
            if nargin < 4
                indices = [];
            end
            
            msmsd = obj.getMeanMSD(indices);
            
            if nargin < 3
                errorbar = false;
                if nargin < 2
                    ha = gca;
                end
            end
            
            if errorbar
                h = msdanalyzer.errorShade(ha, msmsd(:,1), msmsd(:,2), msmsd(:,3), [0 0 0], false);
                set(h.mainLine, 'LineWidth', 2);
                
            else
                h = plot(ha, msmsd(:,1), msmsd(:,2), 'k', ...
                    'LineWidth', 2);
            end
            
            obj.labelPlotMSD(ha);
            
            if nargout > 0
                varargout{1} = h;
                if nargout > 1
                    varargout{2} = ha;
                end
            end
            
        end
        
        
        function obj = computeDrift(obj, method, extra, interpmethod)
            %%COMPUTEDRIFT Compute and store drift correction.
            %
            % obj = obj.computeDrift(method) computes and stores the drift
            % using one of the 4 following methods:
            %
            % 'clear' does not compute drift and remove any prior drift
            % computation results.
            %
            % 'manual' allow to specify manually the drift vector:
            % obj = obj.computeDrift('manual', dv); where dv is a double array
            % of size N x (nDim+1) (nDim being the problem dimensionality), and
            % must be arranged as following: [ Ti Xi Yi ... ] etc...
            % On top of this, the drift vector must cover all the possible time
            % points specified in the tracks field of this object: It must
            % start before the first point and end after the last one,
            % otherwise an error is thrown.
            %
            % Missing values within these extremities are interpolated using a
            % linear interpolation scheme by default. To specify another
            % interpolation scheme, use the following syntax:
            % obj.computeDrift('manual', dv, interpmethod), with interpmethod
            % being any value accepted by interp1 ('linear', 'nearest',
            % 'spline', 'pchip', 'cubic').
            %
            % 'centroid' derives the drift by computing the center of mass of
            % all particles at each time point. This method best work for a
            % large number of particle and when the same number of particles is
            % found at every time points. It fails silently otherwise.
            %
            % 'velocity' derives drift by computing instantaneous velocities
            % and averaging them all together, at each time point. If the
            % particles are in sufficient number, and if their real movement is
            % indeed uncorrelated, the uncorrelated part will cancel when
            % averaging, leaving only the correlated part. We assume this part
            % is due to the drift. This method is more robust than the
            % 'centroid' method against particle disappearance and appearance.
            %
            % Results are stored in the 'drift' field of the returned object.
            % It is a double array of size N x (nDim+1) (nDim being the problem
            % dimensionality), and must be arranged as following: [ Ti Xi Yi ... ]
            % etc. If present, it will by used for any call to computeMSD and
            % computeVCorr methods.
            
            
            % First compute common time points
            time = obj.getCommonTimes();
            n_times = numel(time);
            n_tracks = numel(obj.tracks);
            
            switch lower(method)
                
                case 'manual'
                    
                    if nargin < 4
                        interpmethod = 'linear';
                    end
                    
                    drift_dim = size(extra, 2);
                    if drift_dim ~= obj.n_dim + 1
                        error('msdanalyzer:computeDrift:BadDimensionality', ...
                            'Drift must be of size N x (nDim+1) with [ T0 X0 Y0 ... ] etc...');
                    end
                    
                    uninterpolated_time = extra(:, 1);
                    uninterpolated_drift = extra(:, 2:end);
                    if min(uninterpolated_time) > min(time) || max(uninterpolated_time) < max(time)
                        error('msdanalyzer:computeDrift:BadTimeVector', ...
                            'For manual drift correction, time vector must cover all time vector from all tracks.');
                    end
                    
                    ldrift = interp1(...
                        uninterpolated_time, ...
                        uninterpolated_drift, ...
                        time, interpmethod);
                    
                    obj.drift = [time ldrift];
                    
                case 'centroid'
                    
                    ldrift = zeros(n_times, obj.n_dim);
                    n_drift = zeros(n_times, 1);
                    for i = 1 : n_tracks
                        
                        t = obj.tracks{i}(:,1);
                        t = msdanalyzer.roundn(t, msdanalyzer.TOLERANCE);
                        
                        % Determine target time index in bulk
                        [~, index_in_all_tracks_time, ~] = intersect(time, t);
                        
                        % Add to mean accum for these indexes
                        n_drift(index_in_all_tracks_time) = n_drift(index_in_all_tracks_time) + 1;
                        ldrift(index_in_all_tracks_time, :) = ldrift(index_in_all_tracks_time, :) + obj.tracks{i}(:, 2:end);
                        
                    end
                    
                    ldrift = ldrift ./ repmat(n_drift, [1 obj.n_dim]);
                    obj.drift = [time ldrift];
                    
                    
                case 'velocity'
                    
                    sum_V = zeros(n_times, obj.n_dim);
                    n_V = zeros(n_times, 1);
                    
                    for i = 1 : n_tracks
                        
                        t = obj.tracks{i}(:,1);
                        t = msdanalyzer.roundn(t, msdanalyzer.TOLERANCE);
                        
                        % Determine target time index in bulk
                        [~, index_in_all_tracks_time, ~] = intersect(time, t);
                        
                        % Remove first element
                        index_in_all_tracks_time(1) = [];
                        
                        % Compute speed
                        V = diff( obj.tracks{i}(:, 2:end) ) ./ repmat(diff(t), [ 1 obj.n_dim]);
                        
                        % Add to mean accum for these indexes
                        n_V(index_in_all_tracks_time) = n_V(index_in_all_tracks_time) + 1;
                        sum_V(index_in_all_tracks_time, :) = sum_V(index_in_all_tracks_time, :) + V;
                        
                    end
                    
                    % Build accumulated drift
                    sum_V(1, :) = 0;
                    n_V(1, :) = 1;
                    % Integrate
                    d_time = [0; diff(time) ];
                    ldrift = cumsum( sum_V ./ repmat(n_V, [1 obj.n_dim]) .* repmat(d_time, [1 obj.n_dim]), 1);
                    obj.drift = [time ldrift];
                    
                case 'clear'
                    
                    obj.drift = [];
                    
                    
                otherwise
                    error('msdanalyzer:computeDriftCorrection:UnknownCorrectionMethod', ...
                        'Unknown correction method %s. Must be ''clear'', ''manual'', ''centroid'' or ''velocity''.', ...
                        method);
            end
            
            obj.msd_valid = false;
            obj.vcorr_valid = false;
            
        end
        
        
        function velocities = getVelocities(obj, indices)
            %%GETVELOCITIES Generate and return the instantaneous velocities.
            %
            % v = obj.getVelocities returns in v the instantaneous velocities
            % calculated over all the particles tracjectories stored in this
            % object, using dX/dt. The velocities are corrected for drift if
            % this object holds a proper drift field.
            %
            % This method returns a cell array, one cell per particle. Arrays
            % are N x (Ndim+1) double arrays, with Ndim the dimensionality set
            % at object creation. Data is organized as follow:  [ Ti Vxi Vyi ... ].
            %
            % v = obj.getVelocities(indices) restrict the calculation over only
            % the particles with specified indices. Use an empty array to use
            % take all.
            
            if nargin < 2 || isempty(indices)
                indices = 1 : numel(obj.tracks);
            end
            
            n_tracks = numel(indices);
            velocities = cell(n_tracks, 1);
            
            for i = 1 : n_tracks
                
                index = indices(i);
                
                t = obj.tracks{index}(:, 1);
                X = obj.tracks{index}(:, 2:end);
                
                % Determine drift correction
                if ~isempty(obj.drift)
                    tdrift = obj.drift(:, 1);
                    xdrift = obj.drift(:, 2:end);
                    % Determine target delay index in bulk
                    [~, index_in_drift_time, index_in_track_time] = intersect(tdrift, t);
                    % Keep only track times that can be corrected.
                    X = X(index_in_track_time, :);
                    t = t(index_in_track_time);
                    % Subtract drift position to track position
                    X = X - xdrift(index_in_drift_time, :);
                    
                end
                
                dX = diff(X, 1) ./ repmat(diff(t), [1 obj.n_dim]);
                velocities{i} = [ t(1:end-1) dX];
            end
            
        end
        
        
        function obj = computeVCorr(obj, indices)
            %%COMPUTEVCORR Compute velocity autocorrelation.
            %
            % obj = obj.computeVCorr computes the velocity autocorrelation for all
            % the particles trajectories stored in this object. Velocity
            % autocorrelation is defined as vc(t) = < v(i+t) x v(i) >, the mean
            % being taken over all possible pairs inside a trajectories.
            %
            % Results are stored in the 'vcorr' field of the returned object.
            % The velocity autocorrelation is stored for each particles in a
            % cell array, one cell per particle. The array is a double array of
            % size N x 4, and is arranged as follow: [dt mean std N ; ...]
            % where dt is the delay for the autocorrelation, mean is the mean
            % autocorrelation value for this delay, std the standard deviation
            % and N the number of points in the average.
            %
            % obj = obj.computeVCorr(indices) computes the velocity
            % autocorrelation only for the particles with the specified
            % indices. Use an empty array to take all particles.
            
            obj.vcorr = cell(numel(obj.tracks), 1);
            
            if nargin < 2 || isempty(indices)
                indices = 1 : numel(obj.tracks);
            end
            
            % Get instantaneous velocities
            velocities = obj.getVelocities(indices);
            delays = obj.getAllDelays(indices);
            n_delays = numel(delays);
            n_tracks = numel(velocities);
            
            fprintf('Computing velocity autocorrelation of %d tracks... ', n_tracks);
            fprintf('%4d/%4d', 0, n_tracks);
            for i = 1 : n_tracks
                fprintf('\b\b\b\b\b\b\b\b\b%4d/%4d', i, n_tracks);
                
                % Holder for mean, std calculations
                sum_vcorr     = zeros(n_delays-1, 1);
                sum_vcorr2    = zeros(n_delays-1, 1);
                n_vcorr       = zeros(n_delays-1, 1);
                
                % Unwrap data
                vc = velocities{i};
                t = vc(:, 1);
                V = vc(:, 2:end);
                n_detections = size(V, 1);
                
                % First compute velocity correleation at dt = 0 over all tracks
                Vc0     = mean( sum( V.^2, 2) );
                
                % Other dts
                for j = 1 : n_detections - 1
                    
                    % Delay in physical units
                    dt = t(j+1:end) - t(j);
                    dt = msdanalyzer.roundn(dt, msdanalyzer.TOLERANCE);
                    
                    % Determine target delay index in bulk
                    [~, index_in_all_delays, ~] = intersect(delays, dt);
                    
                    % Velocity correlation in bulk
                    lvcorr = sum( repmat(V(j, :), [ (n_detections-j) 1]) .* V(j+1:end, :), 2 );
                    
                    % Normalize
                    lvcorr = lvcorr ./ Vc0;
                    
                    % Store for mean computation
                    sum_vcorr(index_in_all_delays)   = sum_vcorr(index_in_all_delays) + lvcorr;
                    sum_vcorr2(index_in_all_delays)  = sum_vcorr2(index_in_all_delays) + (lvcorr .^ 2);
                    n_vcorr(index_in_all_delays)     = n_vcorr(index_in_all_delays) + 1;
                    
                end
                
                mean_vcorr = sum_vcorr ./ n_vcorr;
                std_vcorr = sqrt( (sum_vcorr2 ./ n_vcorr) - (mean_vcorr .* mean_vcorr)) ;
                vcorrelation = [ delays(1:end-1) mean_vcorr std_vcorr n_vcorr ];
                vcorrelation(1,:) = [0 1 0 n_detections];
                
                % Store in object field
                index = indices(i);
                obj.vcorr{index} = vcorrelation;
                
            end
            fprintf('\b\b\b\b\b\b\b\b\bDone.\n')
            
            obj.vcorr_valid = true;
            
        end
        
        
        function msvcorr = getMeanVCorr(obj, indices)
            %%GETMEANVCORR Compute the weighted mean of velocity autocorrelation.
            %
            % msd = obj.getMeanVCorr computes and return the weigthed mean of
            % all velocity autocorrelation curves stored in this object. All
            % possible delays are first derived, and for each delay, a weighted
            % mean is computed from all the velocity autocorrelation curves
            % stored in this object. Weights are set to be the number of points
            % averaged to generate the mean square displacement value at the
            % given delay. Thus, we give more weight to velocity
            % autocorrelation curves with greater certainty (larger number of
            % elements averaged).
            %
            % Results are returned as a N x 4 double array, and ordered as
            % following: [ dT M STD N ] with:
            % - dT the delay vector
            % - M the weighted mean of velocity autocorrelation for each delay
            % - STD the weighted standard deviation
            % - N the number of degrees of freedom in the weighted mean
            % (see http://en.wikipedia.org/wiki/Weighted_mean)
            %
            % msd = obj.getMeanVCorr(indices) only takes into account the
            % velocity autocorrelation curves with the specified indices,
            
            if ~obj.vcorr_valid
                obj = obj.computeVCorr(indices);
            end
            
            if nargin < 2 || isempty(indices)
                indices = 1 : numel(obj.vcorr);
            end
            
            n_tracks = numel(indices);
            
            % First, collect all possible delays
            all_delays = cell(n_tracks, 1);
            for i = 1 : n_tracks
                index = indices(i);
                all_delays{i} = obj.vcorr{index}(:,1);
            end
            delays = unique( vertcat( all_delays{:} ) );
            n_delays = numel(delays);
            
            
            % Collect
            sum_weight          = zeros(n_delays, 1);
            sum_weighted_mean   = zeros(n_delays, 1);
            
            % 1st pass
            for i = 1 : n_tracks
                
                index = indices(i);
                
                t = obj.vcorr{index}(:,1);
                m = obj.vcorr{index}(:,2);
                n = obj.vcorr{index}(:,4);
                
                % Do not tak NaNs
                valid = ~isnan(m);
                t = t(valid);
                m = m(valid);
                n = n(valid);
                
                % Find common indices
                [~, index_in_all_delays, ~] = intersect(delays, t);
                
                % Accumulate
                sum_weight(index_in_all_delays)           = sum_weight(index_in_all_delays)         + n;
                sum_weighted_mean(index_in_all_delays)    = sum_weighted_mean(index_in_all_delays)  + m .* n;
            end
            
            % Compute weighted mean
            mmean = sum_weighted_mean ./ sum_weight;
            
            % 2nd pass: unbiased variance estimator
            sum_weighted_variance = zeros(n_delays, 1);
            sum_square_weight     = zeros(n_delays, 1);
            
            for i = 1 : n_tracks
                
                index = indices(i);
                
                t = obj.vcorr{index}(:,1);
                m = obj.vcorr{index}(:,2);
                n = obj.vcorr{index}(:,4);
                
                % Do not tak NaNs
                valid = ~isnan(m);
                t = t(valid);
                m = m(valid);
                n = n(valid);
                
                % Find common indices
                [~, index_in_all_delays, ~] = intersect(delays, t);
                
                % Accumulate
                sum_weighted_variance(index_in_all_delays)    = sum_weighted_variance(index_in_all_delays)  + n .* (m - mmean(index_in_all_delays)).^2 ;
                sum_square_weight(index_in_all_delays)        = sum_square_weight(index_in_all_delays)      + n.^2;
            end
            
            % Standard deviation
            mstd = sqrt( sum_weight ./ (sum_weight.^2 - sum_square_weight) .* sum_weighted_variance );
            
            % Output [ T mean std Nfreedom ]
            msvcorr = [ delays mmean mstd (sum_weight.^2 ./ sum_square_weight) ];
            
        end
        
        
        function varargout = plotMeanVCorr(obj, ha, errorbar, indices)
            %%PLOTMEANVCORR Plot the weighted mean of the velocity autocorrelation curves.
            %
            % obj,plotMeanVCorr computes and plots the weighted of all velocity
            % autocorrelation curves. See msdanalyzer.getMeanVCorr.
            %
            % obj,plotMeanVCorr(ha) plots the curve in the axes with the
            % specified handle.
            %
            % obj,plotMeanVCorr(ha, errorbar) where 'errorbar' is a boolean
            % allow to specify whether to plot the curve with error bars
            % indicating the weighted standard deviation. Default is false.
            %
            % obj,plotMeanVCorr(ha, errorbar, indices) computes and plots the
            % mean only fothe velocity autocorrelation curves whose indices are
            % given on the 'indices' array.
            %
            % h = obj,plotMeanVCorr(...) returns the handle to the line plotted.
            %
            % [h, ha] = obj,plotMeanVCorr(...) also returns the handle of the
            % axes in which the curve was plotted.
            
            if nargin < 4
                indices = [];
                
                if nargin < 3
                    errorbar = false;
                    if nargin < 2
                        ha = gca;
                    end
                end
            end
            mvc = obj.getMeanVCorr(indices);
            
            if errorbar
                h = msdanalyzer.errorShade(ha, mvc(:,1), mvc(:,2), mvc(:,3), [0 0 0], false);
                set(h.mainLine, 'LineWidth', 2);
                
            else
                h = plot(ha, mvc(:,1), mvc(:,2), 'k', ...
                    'LineWidth', 2);
            end
            
            obj.labelPlotVCorr(ha);
            
            if nargout > 0
                varargout{1} = h;
                if nargout > 1
                    varargout{2} = ha;
                end
            end
            
        end
        
    end
    
    %% Private methods
    
    methods (Access = private)
        
        function time = getCommonTimes(obj)
            
            n_tracks = numel(obj.tracks);
            times = cell(n_tracks, 1);
            for i = 1 : n_tracks
                times{i} = obj.tracks{i}(:,1);
            end
            time = unique( vertcat(times{:}) );
            time = msdanalyzer.roundn(time, msdanalyzer.TOLERANCE);
            
        end
        
        function delays = getAllDelays(obj, indices)
            % First, find all possible delays in time vectors.
            % Time can be arbitrary spaced, with frames missings,
            % non-uniform sampling, etc... so we have to do this clean.
            
            if nargin < 2 || isempty(indices)
                indices = 1 : numel(obj.tracks);
            end
            
            n_tracks = numel(indices);
            all_delays = cell(n_tracks, 1);
            for i = 1 : n_tracks
                index = indices(i);
                track = obj.tracks{index};
                t = track(:,1);
                [T1, T2] = meshgrid(t, t);
                dT = msdanalyzer.roundn(abs(T1(:)-T2(:)), msdanalyzer.TOLERANCE);
                all_delays{i} = unique(dT);
            end
            delays = unique( vertcat(all_delays{:}) );
        end
        
    end
    
    %% Static methods
    
    methods (Static, Access = private)
        
        function wm = weightedmean(x, w)
            wm = sum( x .* w) ./ sum(w);
        end
        
        function sewm = standarderrorweightedmean(x, w)
            n = numel(w);
            wbar = mean(w);
            xbar = sum( x .* w) ./ sum(w);
            sewm = n /((n-1) * sum(w)^2) * (sum( (w.*x - wbar*xbar).^2) ...
                - 2 * xbar * sum( (w-wbar).*(w.*x - wbar*xbar)) ...
                + xbar^2 * sum((w-wbar).^2));
        end
        
        
        
        function H = errorShade(ha, x, y, errBar, col, transparent)
            % Adapted from Rob Campbell code, at:
            % http://www.mathworks.com/matlabcentral/fileexchange/26311-shadederrorbar/content/shadedErrorBar.m
            hold on
            H.mainLine = plot(ha, x, y, 'Color', col);
            
            edgeColor = col + (1-col) * 0.55;
            patchSaturation = 0.15; %How de-saturated or transparent to make the patch
            if transparent
                faceAlpha=patchSaturation;
                patchColor=col;
                set(gcf,'renderer','openGL')
            else
                faceAlpha=1;
                patchColor=col+(1-col)*(1-patchSaturation);
                set(gcf,'renderer','painters')
            end
            
            %Calculate the y values at which we will place the error bars
            uE = y + errBar;
            lE = y - errBar;
            
            %Make the cordinats for the patch
            yP = [ lE ; flipud(uE) ];
            xP = [ x ; flipud(x) ];
            
            invalid = isnan(xP) | isnan(yP) | isinf(xP) | isinf(yP);
            yP(invalid) = [];
            xP(invalid) = [];
            
            
            H.patch = patch(xP, yP, 1, ...
                'Facecolor', patchColor,...
                'Edgecolor', 'none',...
                'Facealpha', faceAlpha, ...
                'Parent', ha);
            
            %Make nice edges around the patch.
            H.edge(1) = plot(ha, x, lE, '-', 'Color', edgeColor);
            H.edge(2) = plot(ha, x, uE, '-', 'Color', edgeColor);
            
            %The main line is now covered by the patch object and was plotted first to
            %extract the RGB value of the main plot line. I am not aware of an easy way
            %to change the order of plot elements on the graph so we'll just remove it
            %and put it back (yuk!)
            delete(H.mainLine)
            H.mainLine = plot(ha, x, y, 'Color', col);
        end
        
    end
    
    %% Static and public functions
    
    methods (Static)
        
        %ROUNDN  Round towards nearest number with Nth decimals.
        %   ROUNDN(X,N) rounds the elements of X to the nearest numbers with the
        %   precision given by N.
        %
        %   Examples:   roundn(8.73,0) = 9
        %               roundn(8.73,1) = 8.7
        %               roundn(8.73,2) = 8.73
        %               roundn(8.73,-1) = 10
        %
        %   See also ROUND
        % Jean-Yves Tinevez - MPI-CBG - August 2007
        
        function Y = roundn(X,N)
            Y = 10^(-N) .* round(X.*10^(N));
        end
        
        function newobj = pool(msda_arr)
            %%POOL Pool the data of several masanalyzer objects in a new one.
            
            if ~isa(msda_arr, 'msdanalyzer')
                error('msdanalyzer:pool:BadArgument', ...
                    'Expect arguments to be of class ''masanalyzer'', got ''%s''.',...
                    class(msda_arr))
            end
            
            n_obj = numel(msda_arr);
            
            % Check dimensionality and units consistency
            lspace_units = msda_arr(1).space_units;
            ltime_units = msda_arr(1).time_units;
            ln_dim = msda_arr(1).n_dim;
            
            for i = 2 : n_obj
                
                obj = msda_arr(i);
                if ~strcmp(obj.space_units, lspace_units)
                    error('msdanalyzer:pool:InconsistentArray', ...
                        'Element %d is inconsistent. Expected space units to be %s, got %s.', ...
                        i, lspace_units, obj.space_units)
                end
                if ~strcmp(obj.time_units, ltime_units)
                    error('msdanalyzer:pool:InconsistentArray', ...
                        'Element %d is inconsistent. Expected time units to be %s, got %s.', ...
                        i, ltime_units, obj.time_units)
                end
                if obj.n_dim ~= ln_dim
                    error('msdanalyzer:pool:InconsistentArray', ...
                        'Element %d is inconsistent. Expected dimensionality to be %d, got %d.', ...
                        i, ln_dim, obj.n_dim)
                end
                
            end
            
            all_msd = cell(n_obj,1 );
            all_vcorr = cell(n_obj,1 );
            
            for i = 1 : n_obj
                obj = msda_arr(i);
                
                obj = obj.computeDrift('velocity');
                obj = obj.computeMSD;
                obj = obj.computeVCorr;
                
                all_msd{i} = obj.msd;
                all_vcorr{i} = obj.vcorr;
            end
            
            newobj = msdanalyzer(ln_dim, lspace_units, ltime_units);
            newobj.msd = vertcat(all_msd{:});
            newobj.vcorr = vertcat(all_vcorr{:});
            newobj.msd_valid = true;
            newobj.vcorr_valid = true;
            
        end
        
    end
    
end

Contact us