function [pivDataOut] = pivManipulateData(action,varargin)
% pivManipulateData - manipulation of data structures generated by PIVsuite
%
% Usage:
%    [pivDataOut] = pivManipulateData(action,varargin)
%
%    Typical usage is:
%    1. pivDataSeq = pivManipulateData('initSequenceData',pivData,nt)
%          Input pivData is a result of processing of a single pair of images (fields pivData.U, .V, .Status
%          etc. are 2D arrays). Output pivDataSeq then contains these arrays with added third dimension for
%          time (with size given by nt, which is number of image pairs). Input pivData is copied in this
%          structure as a result of the first image pair.
%        - Note: fields pivData.spuriousU, .spuriousV, .spuriousX, .spuriousY are not copied. X and Y
%          information about spurious vectors can be obtained from the .Status field.
%    2. pivDataSeq = pivManipulateData('writeTimeSlice',pivDataSeq,pivData,n)
%          Input pivDataSeq is is the data structure for a sequence of images (see usage 1 for its
%          initialization). Input pivData contains result for n-th pair of images (/n-th time slice).
%          pivManipulateData inserts pivData to n-th position of pivDataSeq.
%    3. pivData = pivManipulateData('readTimeSlice',pivDataSeq,n)
%          returns n-th time slice contained in pivDataSeq
%    4. pivData = pivManipulateData('multiplyVelocity', pivData, mult)
%          Returns pivData, in which all velocity information is multiplied by factor mult.
%    5. pivDataSeq = pivManipulateData('multiplyVelocity', pivDataSeq, mult)
%          Returns pivDataSeq, in which all velocity information is multiplied by factor mult.
%    6. pivData = pivManipulateData('combineData',pivData1,pivData2)
%       pivData = pivManipulateData('combineData',pivData1,pivData2,UmagMax)
%       pivData = pivManipulateData('combineData',pivData1,pivData2,UmagMax,minDetectability)
%       pivData = pivManipulateData('combineData',pivData1,pivData2,UmagMax,minDetectability,mult1,mult2)
%          Combines to velocity fields measured with different time-step. It is assumed that time-step of
%          pivDataX is multX times longer than the interval between subsequent frames. The pivDataX with
%          smaller multX factor is considered as less precise (smaller displacement), but more robust (less
%          spurious vectors). Following rules are adopted:
%             - take the data of velocity field with higher mult factor as default output
%             - if .Status of data with higher multX indicates either cross-correlation failure or spurious
%               vector, replace them with velocity with smaller mult
%             - if quality score of data with smaller multX factor is better than score of those with higher
%               mult factor, replace them with velocity with smaller mult
%             - quality score of each velocity vector is
%                  score = mult * min(Umag,UmagMax) * max((ccPeak-ccPeaksecondary-minDetectability), 0).
%               If Umag is larger than UmagMax, it is clipped to this value.
%        - If maximum velocity magnitude is not specified, it is set to default value UmagMax = (iaSizeX +
%          iaSizeY) / 2. If minDetectability is not specified, its default value is 0.08.
%        - Note: fields pivData.spuriousU, .spuriousV, .spuriousX, .spuriousY are not processed.
%        - Note: Status is updated from the corresponding velocity field.
%        - Additional field, pivData.multiplicator, is added. This field contains multX factor for each
%          velocity vector in the output field. This field can be used to identify, if the velocity vector
%          comes from pivData1 or pivData2.
%    7. pivDataSeq = pivManipulateData('combineData',pivDataSeq1,pivDataSeq2,mult1,mult2)
%       pivDataSeq = pivManipulateData('combineData',pivDataSeq1,pivDataSeq2,mult1,mult2,UmagMax)
%          Same as Usage 6., but works with pivDataSeq instead of pivData
%    8. pivData = pivManipulateData('limitX',pivData,[xmin,xmax])
%       pivDataSeq = pivManipulateData('limitX',pivDataSeq,[xmin,xmax])
%       pivData = pivManipulateData('limitY',pivData,[ymin,ymax])
%       pivDataSeq = pivManipulateData('limitY',pivDataSeq,[ymin,ymax])
%       pivData = pivManipulateData('limitXY',pivData,[ymin,ymax])
%       pivDataSeq = pivManipulateData('limitXY',pivDataSeq,[xmin,xmax,ymin,ymax])
%          Removes from data all data points, which have X coordinates outside range <xmin,xmax> and/or which
%          have Y coordinates outside range <ymin,ymax>.
%    9. pivDataStat = pivManipulateData('writeToStat',pivDataStat,pivData)
%          Takes the velocity field contained in pivData, and adds it to velocity statistics included in
%          pivDataStat. #################### BETTER DESCRIPTION #######################
%   10. pivDataStat = pivManipulateData('initStat',pivDataStat,pivData)
%          Initializes the data structure for, which will contain velocity statistics.
%          #################### BETTER DESCRIPTION #######################
%
% Inputs:
%    action ... defines, how the data will be manipulated. Recognized actions (case-insensitive) are
%         'initSequenceData', 'writeTimeSlice', 'readTimeSlice','multiplyVelocity','tSingle','toDouble'.
%         See "Usage" for details.
%    pivData ... structure containing results of pivAnalyzeImagePair.m. In this structure, velocity fields .U
%         and .V have ny x nx elements (where ny and nx is the number of interrogation areas).
%    pivDataSeq ... structure containing results of pivAnalyzeImageSequence.m. In this structure, velocity
%         fields .U and .V have ny x nx x nt elements (where ny and nx is the number of interrogation areas,
%         nt is the number of image pairs). .ccPeak and .ccPeakSecondary have the same size. Fields
%         .spuriousU, .spuriousV, .spuriousX, .spuriousY are missing.
%    n ... number of time slice
%    nt ... total number of time slices
%    mult ... velocity multiplier
%    UmagMax ... maximum velocity amplitude; above this level, quality score is not further increased when
%         velocity magnitude is increasing.
%
% Outputs:
%    pivData, pivDataSeq ... see inputs for their meaning.
%
%
%
% This subroutine is a part of
%
% =========================================
%               PIVsuite
% =========================================
%
% PIVsuite is a set of subroutines intended for processing of data acquired with PIV (particle image
% velocimetry).
%
% Written by Jiri Vejrazka, Institute of Chemical Process Fundamentals, Prague, Czech Republic
%
% For the use, see files example_XX_xxxxxx.m, which acompany this file. PIVsuite was tested with
% Matlab 7.12 (R2011a) and 7.14 (R2012a).
%
% In the case of a bug, please, contact me: vejrazka (at) icpf (dot) cas (dot) cz
%
%
% Requirements:
%     Image Processing Toolbox
%
%     inpaint_nans.m
%         subroutine by John D'Errico, available at http://www.mathworks.com/matlabcentral/fileexchange/4551
%
%     smoothn.m
%         subroutine by Damien Garcia, available at
%         http://www.mathworks.com/matlabcentral/fileexchange/274-smooth
%
% Credits:
%    PIVsuite is a redesigned version of PIVlab software [3], developped by W. Thielicke and E. J. Stamhuis.
%    Some parts of this code are copied or adapted from it (especially from its piv_FFTmulti.m subroutine).
%    PIVsuite uses 3rd party software:
%        inpaint_nans.m, by J. D'Errico, [2]
%        smoothn.m, by Damien Garcia, [5]
%
% References:
%   [1] Adrian & Whesterweel, Particle Image Velocimetry, Cambridge University Press 2011
%   [2] John D'Errico, inpaint_nans subroutine, http://www.mathworks.com/matlabcentral/fileexchange/4551
%   [3] W. Thielicke and E. J. Stamhuid, PIVlab 1.31, http://pivlab.blogspot.com
%   [4] Raffel, Willert, Wereley & Kompenhans, Particle Image Velocimetry: A Practical Guide. 2nd edition,
%       Springer 2007
%   [5] Damien Garcia, smoothn subroutine, http://www.mathworks.com/matlabcentral/fileexchange/274-smooth
%
%
% Acronyms and meaning of variables used in this subroutine:
%    IA ... concerns "Interrogation Area"
%    im ... image
%    dx ... some index
%    ex ... expanded (image)
%    est ... estimate (velocity from previous pass) - will be used to deform image
%    aux ... auxiliary variable (which is of no use just a few lines below)
%    cc ... cross-correlation
%    vl ... validation
%    sm ... smoothing
%    Word "velocity" should be understood as "displacement"



switch(lower(action))
    % ACTION initSequenceData
    case 'initsequencedata'
        pivDataIn = varargin{1};
        nt = varargin{2};           % number of time slices
        % get the size of measurement field
        nx = size(pivDataIn.X,2);
        ny = size(pivDataIn.X,1);
        % copy information to output structure
        pivDataOut.X = single(pivDataIn.X);
        pivDataOut.Y = single(pivDataIn.Y);
        pivDataOut.U = single(zeros(ny,nx,nt) + NaN);
        pivDataOut.U(:,:,1) = single(pivDataIn.U);
        pivDataOut.V = single(zeros(ny,nx,nt) + NaN);
        pivDataOut.V(:,:,1) = single(pivDataIn.V);
        pivDataOut.Status = uint16(zeros(ny,nx,nt)+255);
        pivDataOut.Status(:,:,1) = uint16(pivDataIn.Status);
        pivDataOut.N = pivDataIn.N;
        pivDataOut.iaSizeX = pivDataIn.iaSizeX;
        pivDataOut.iaSizeY = pivDataIn.iaSizeY;
        pivDataOut.iaStepX = pivDataIn.iaStepX;
        pivDataOut.iaStepY = pivDataIn.iaStepY;
        pivDataOut.imSizeX = pivDataIn.imSizeX;
        pivDataOut.imSizeY = pivDataIn.imSizeY;
        pivDataOut.ccPeak = single(zeros(ny,nx,nt) + NaN);
        pivDataOut.ccPeak(:,:,1) = single(pivDataIn.ccPeak);
        pivDataOut.ccPeakSecondary = single(zeros(ny,nx,nt) + NaN);
        pivDataOut.ccPeakSecondary(:,:,1) = single(pivDataIn.ccPeakSecondary);
        if isfield(pivDataIn,'ccStd1')
            pivDataOut.ccStd1 = single(zeros(ny,nx,nt) + NaN);
            pivDataOut.ccStd1(:,:,1) = single(pivDataIn.ccStd1);
        end
        if isfield(pivDataIn,'ccStd2')
            pivDataOut.ccStd2 = single(zeros(ny,nx,nt) + NaN);
            pivDataOut.ccStd2(:,:,1) = single(pivDataIn.ccStd2);
        end
        if isfield(pivDataIn,'ccMean1')
            pivDataOut.ccMean1 = single(zeros(ny,nx,nt) + NaN);
            pivDataOut.ccMean1(:,:,1) = single(pivDataIn.ccMean1);
        end
        if isfield(pivDataIn,'ccMean2')
            pivDataOut.ccMean2 = single(zeros(ny,nx,nt) + NaN);
            pivDataOut.ccMean2(:,:,1) = single(pivDataIn.ccMean2);
        end
        pivDataOut.ccFailedN = zeros(nt,1);
        pivDataOut.ccFailedN(1) = pivDataIn.ccFailedN;
        pivDataOut.ccSubpxFailedN = zeros(nt,1);
        pivDataOut.ccSubpxFailedN(1) = pivDataIn.ccSubpxFailedN;
        pivDataOut.maskedN = zeros(nt,1);
        pivDataOut.maskedN(1) = pivDataIn.maskedN;
        pivDataOut.replacedN = zeros(nt,1);
        pivDataOut.replacedN(1) = pivDataIn.replacedN;
        pivDataOut.spuriousN = zeros(nt,1);
        pivDataOut.spuriousN(1) = pivDataIn.spuriousN;
        if isfield(pivDataIn,'imFilename1')
            pivDataOut.imFilename1 = cell(nt,1);
            pivDataOut.imFilename1{1,1} = pivDataIn.imFilename1;
        end
        if isfield(pivDataIn,'imMaskFilename1')
            pivDataOut.imMaskFilename1 = cell(nt,1);
            pivDataOut.imMaskFilename1{1,1} = pivDataIn.imMaskFilename1;
        end
        if isfield(pivDataIn,'imNo1')
            pivDataOut.imNo1 = zeros(nt,1)+NaN;
            pivDataOut.imNo1(1) = pivDataIn.imNo1;
        end
        if isfield(pivDataIn,'imFilename2')
            pivDataOut.imFilename2 = cell(nt,1);
            pivDataOut.imFilename2{1,1} = pivDataIn.imFilename2;
        end
        if isfield(pivDataIn,'imMaskFilename2')
            pivDataOut.imMaskFilename2 = cell(nt,1);
            pivDataOut.imMaskFilename2{1,1} = pivDataIn.imMaskFilename2;
        end
        if isfield(pivDataIn,'imNo2')
            pivDataOut.imNo2 = zeros(nt,1)+NaN;
            pivDataOut.imNo2(1) = pivDataIn.imNo2;
        end
        if isfield(pivDataIn,'imPairNo')
            pivDataOut.imPairNo = zeros(nt,1)+NaN;
            pivDataOut.imPairNo(1) = pivDataIn.imPairNo;
        end
        if isfield(pivDataIn,'multiplicator')
            pivDataOut.multiplicator = zeros(ny,nx,nt)+NaN;
            pivDataOut.multiplicator(:,:,1) = pivDataIn.multiplicator;
        end
        
        
    % ACTION writeTimeSlice
    case 'writetimeslice'
        pivDataOut = varargin{1};
        pivDataIn = varargin{2};
        kt = varargin{3};         %timeslice number
        % Check, if X and Y are same. If not, give a message and quit.
        if ~(isequal(pivDataOut.X,pivDataIn.X)&&isequal(pivDataOut.Y,pivDataIn.Y))
            disp('Error (pivManipulate.m, action "writeTimeSlice"): Measurement grid is not compatible. Quitting.');
            return;
        end
        % copy information to output structure
        pivDataOut.U(:,:,kt) = single(pivDataIn.U);
        pivDataOut.V(:,:,kt) = single(pivDataIn.V);
        pivDataOut.Status(:,:,kt) = uint16(pivDataIn.Status);
        pivDataOut.ccPeak(:,:,kt) = single(pivDataIn.ccPeak);
        pivDataOut.ccPeakSecondary(:,:,kt) = single(pivDataIn.ccPeakSecondary);
        if isfield(pivDataIn,'ccStd1')
            pivDataOut.ccStd1(:,:,kt) = single(pivDataIn.ccStd1);
        end
        if isfield(pivDataIn,'ccStd2')
            pivDataOut.ccStd2(:,:,kt) = single(pivDataIn.ccStd2);
        end
        if isfield(pivDataIn,'ccMean1')
            pivDataOut.ccMean1(:,:,kt) = single(pivDataIn.ccMean1);
        end
        if isfield(pivDataIn,'ccMean2')
            pivDataOut.ccMean2(:,:,kt) = single(pivDataIn.ccMean2);
        end
        pivDataOut.N = pivDataIn.N;
        pivDataOut.iaSizeX = pivDataIn.iaSizeX;
        pivDataOut.iaSizeY = pivDataIn.iaSizeY;
        pivDataOut.iaStepX = pivDataIn.iaStepX;
        pivDataOut.iaStepY = pivDataIn.iaStepY;
        pivDataOut.imSizeX = pivDataIn.imSizeX;
        pivDataOut.imSizeY = pivDataIn.imSizeY;
        pivDataOut.ccFailedN(kt) = pivDataIn.ccFailedN;
        pivDataOut.ccSubpxFailedN(kt) = pivDataIn.ccSubpxFailedN;
        pivDataOut.maskedN(kt) = pivDataIn.maskedN;
        pivDataOut.replacedN(kt) = pivDataIn.replacedN;
        pivDataOut.spuriousN(kt) = pivDataIn.spuriousN;
        if isfield(pivDataIn,'imFilename1')
            pivDataOut.imFilename1{kt,1} = pivDataIn.imFilename1;
        end
        if isfield(pivDataIn,'imMaskFilename1')
            pivDataOut.imMaskFilename1{kt,1} = pivDataIn.imMaskFilename1;
        end
        if isfield(pivDataIn,'imNo1')
            pivDataOut.imNo1(kt) = pivDataIn.imNo1;
        end
        if isfield(pivDataIn,'imFilename2')
            pivDataOut.imFilename2{kt,1} = pivDataIn.imFilename2;
        end
        if isfield(pivDataIn,'imMaskFilename2')
            pivDataOut.imMaskFilename2{kt,1} = pivDataIn.imMaskFilename2;
        end
        if isfield(pivDataIn,'imNo2')
            pivDataOut.imNo2(kt) = pivDataIn.imNo2;
        end
        if isfield(pivDataIn,'imPairNo')
            pivDataOut.imPairNo(kt) = pivDataIn.imPairNo;
        end
        if isfield(pivDataIn,'multiplicator')
            pivDataOut.multiplicator(:,:,kt) = pivDataIn.multiplicator;
        end
        % copy all fields starting with 'exp'
        fields = fieldnames(pivDataIn);
        auxI = strfind(fields,'exp');
        for kk = 1:numel(auxI)
            if numel(auxI{kk})>0 && auxI{kk,1}==1
                pivDataOut.(fields{kk}) = pivDataIn.(fields{kk});
            end
        end
        
        
    % ACTION readTimeSlice
    case 'readtimeslice'
        pivDataIn = varargin{1};
        kt = varargin{2};         %timeslice number
        % copy information to output structure
        pivDataOut = pivDataIn;
        pivDataOut.U = double(pivDataIn.U(:,:,kt));
        pivDataOut.V = double(pivDataIn.V(:,:,kt));
        pivDataOut.Status = uint16(pivDataIn.Status(:,:,kt));
        pivDataOut.ccPeak = double(pivDataIn.ccPeak(:,:,kt));
        pivDataOut.ccPeakSecondary = double(pivDataIn.ccPeakSecondary(:,:,kt));
        if isfield(pivDataOut,'vorticity')
            pivDataOut.vorticity = single(pivDataIn.vorticity(:,:,kt));
        end
        if isfield(pivDataOut,'epsLEPIV')
            pivDataOut.epsLEPIV = single(pivDataIn.epsLEPIV(:,:,kt));
        end
        if isfield(pivDataOut,'ccStd1')
            pivDataOut.ccStd1 = single(pivDataIn.ccStd1(:,:,kt));
        end
        if isfield(pivDataOut,'ccStd2')
            pivDataOut.ccStd2 = single(pivDataIn.ccStd2(:,:,kt));
        end
        if isfield(pivDataOut,'ccMean1')
            pivDataOut.ccMean1 = single(pivDataIn.ccMean1(:,:,kt));
        end
        if isfield(pivDataOut,'ccMean2')
            pivDataOut.ccMean2 = single(pivDataIn.ccMean2(:,:,kt));
        end
        if isfield(pivDataOut,'imFilename1')
            if iscell(pivDataOut.imFilename1)
                pivDataOut.imFilename1 = pivDataOut.imFilename1{kt};
            end
        end
        if isfield(pivDataOut,'imFilename2')
            if iscell(pivDataOut.imFilename2)
                pivDataOut.imFilename2 = pivDataOut.imFilename2{kt};
            end
        end
        
    % ACTION initStatData
    case 'initstatdata'
        pivDataOut = varargin{1};
        pivDataIn = varargin{2};
        % get the size of measurement field
        nx = size(pivDataIn.X,2);
        ny = size(pivDataIn.X,1);
        % copy information to output structure
        pivDataOut.X = single(pivDataIn.X);
        pivDataOut.Y = single(pivDataIn.Y);
        pivDataOut.statU = zeros(ny,nx);
        pivDataOut.statU2 = zeros(ny,nx);
        pivDataOut.statV = zeros(ny,nx);
        pivDataOut.statV2 = zeros(ny,nx);
        pivDataOut.statUV = zeros(ny,nx);
        pivDataOut.statN = zeros(ny,nx);
        pivDataOut.Umean = zeros(ny,nx);
        pivDataOut.Vmean = zeros(ny,nx);
        pivDataOut.RSuu = zeros(ny,nx);
        pivDataOut.RSvv = zeros(ny,nx);
        pivDataOut.RSuv = zeros(ny,nx);
        pivDataOut.Kuv = zeros(ny,nx);
        pivDataOut.Kuuv = zeros(ny,nx);
        pivDataOut.Kuvv = zeros(ny,nx);
        pivDataOut.iaSizeX = pivDataIn.iaSizeX;
        pivDataOut.iaSizeY = pivDataIn.iaSizeY;
        pivDataOut.iaStepX = pivDataIn.iaStepX;
        pivDataOut.iaStepY = pivDataIn.iaStepY;
        pivDataOut.imSizeX = pivDataIn.imSizeX;
        pivDataOut.imSizeY = pivDataIn.imSizeY;
        
        
    % ACTION writeToStat
    case 'writetostat'
        pivDataOut = varargin{1};
        pivDataIn = varargin{2};

        % reset velocity vectors with untrusted data
        auxNOK = (bitget(pivDataIn.Status,1)+bitget(pivDataIn.Status,2)+bitget(pivDataIn.Status,3)+...
            bitget(pivDataIn.Status,4)+bitget(pivDataIn.Status,7))>0;
             % includes "masked", "cc failed", "peak-detection failed", "spurious" and "spurious in sequence"
             % flags
        auxOK = ~auxNOK;
        pivDataIn.U(auxNOK) = 0;
        pivDataIn.V(auxNOK) = 0;
        % add velocity information to velocity statistics
        pivDataOut.statU = pivDataOut.statU+pivDataIn.U;
        pivDataOut.statU2 = pivDataOut.statU2+pivDataIn.U.^2;
        pivDataOut.statV = pivDataOut.statU+pivDataIn.V;
        pivDataOut.statV2 = pivDataOut.statU2+pivDataIn.V.^2;
        pivDataOut.statUV = pivDataOut.statUV+pivDataIn.U.*pivDataIn.V;
        pivDataOut.statN = pivDataOut.statN + auxOK;
        % compute mean velocity from the statistics
        pivDataOut.Umean = pivDataOut.statU ./ pivDataOut.statN;
        pivDataOut.Vmean = pivDataOut.statV ./ pivDataOut.statN;        
        pivDataOut.RSuu = pivDataOut.statU2 ./ pivDataOut.statN - pivDataOut.Umean.^2;
        pivDataOut.RSvv = pivDataOut.statV2 ./ pivDataOut.statN - pivDataOut.Vmean.^2;
        pivDataOut.RSuv = pivDataOut.statUV ./ pivDataOut.statN - pivDataOut.Umean.*pivDataOut.Vmean;
        pivDataOut.Kuv = 1/2 * (3/2*pivDataOut.RSuu + 3/2*pivDataOut.RSvv);
        pivDataOut.Kuuv = 1/2 * (2*pivDataOut.RSuu + pivDataOut.RSvv);
        pivDataOut.Kuvv = 1/2 * (pivDataOut.RSuu + 2*pivDataOut.RSvv);
        
        
    % ACTION LEPIVdissip
    case 'lepivdissip'
        pivDataOut = varargin{1};
        pivDataIn = varargin{1};
        if nargin>2  % Smagorinsky constant can be provided as an input
            Cs = varargin{2};   % Smagorinsky constant
        else
            Cs = 0.17; 
        end    
        % compute derivatives of the velocity field
        dUdX = diffField(pivDataIn.U,pivDataIn.X,2);
        dUdY = diffField(pivDataIn.U,pivDataIn.Y,1);
        dVdX = diffField(pivDataIn.V,pivDataIn.X,2);
        dVdY = diffField(pivDataIn.V,pivDataIn.Y,1);
        % compute strain rate (eq. 18 in de Jong, Exp Fluids 2009)
        SijSij = 2*dUdX.^2 + 2*dVdY.^2 + 2*dUdX.*dVdY + 3/2*(dUdY+dVdX).^2;
        % compute dissiption rate (eq. 15, 17 and 18 combined, de Jong, Exp Fluids 2009)
        epsilon = 2*Cs^2 * pivDataIn.iaSizeX*pivDataIn.iaSizeY * sqrt(2*SijSij) .* SijSij;
        % compute mean dissipation rate
        epsilonMean = mean(epsilon,3);
        % output results
        pivDataOut.dUdX = dUdX;
        pivDataOut.dUdY = dUdY;
        pivDataOut.dVdX = dVdX;
        pivDataOut.dVdY = dVdY;
        pivDataOut.SijSij = SijSij;
        pivDataOut.epsLEPIV = epsilon;
        pivDataOut.epsMeanLEPIV = epsilonMean;
        
        

        
        
        % ACTION multiplyVelocity
    case 'multiplyvelocity'
        pivDataIn = varargin{1};
        mult = varargin{2};
        % copy everything
        pivDataOut = pivDataIn;
        % multiply velocities
        pivDataOut.U = pivDataOut.U * mult;
        pivDataOut.V = pivDataOut.V * mult;
        if isfield(pivDataIn,'spuriousU')
            pivDataOut.spuriousU = pivDataOut.spuriousU * mult;
        end
        if isfield(pivDataIn,'spuriousV')
            pivDataOut.spuriousV = pivDataOut.spuriousV * mult;
        end
        
        
    % ACTION combineData
    case 'combinedata'
        % assign inputs
        pivData1 = varargin{1};
        pivData2 = varargin{2};
        % Check, if X and Y are same. If not, give a message and quit.
        if ~(isequal(pivData1.X,pivData2.X)&&isequal(pivData1.Y,pivData2.Y))
            disp('Error (pivManipulate.m, action "combineData"): Measurement grid is not compatible. Quitting.');
            return;
        end
        % read optional parameters
        if nargin > 3, UmagMax = varargin{3};
        else UmagMax = (pivData1.iaSizeX + pivData1.iaSizeY + pivData2.iaSizeX + pivData2.iaSizeY)/4;
        end
        if nargin > 4, minDetecteability = varargin{4}; else minDetecteability = 0.08; end
        if nargin > 5, mult1 = varargin{5};
        else
            if isfield(pivData1,'imNo1') && isfield(pivData1,'imNo2')
                mult1 = pivData1.imNo2(1) - pivData1.imNo1(1);
            else
                disp('Error (pivManipulate.m, action "combineData"): Unknown velocity multiplicator for pivData1');
                return
            end
        end
        if nargin > 6, mult2 = varargin{6};
        else
            if isfield(pivData2,'imNo1') && isfield(pivData2,'imNo2')
                mult2 = pivData2.imNo2(1) - pivData2.imNo1(1);
            else
                disp('Error (pivManipulate.m, action "combineData"): Unknown velocity multiplicator for pivData2');
                return
            end
        end
        % Check, for which times the time-slices in pivData1 and pivData2 are. If not for same times, skip them
        auxI1 = (1:numel(pivData1.imPairNo))';
        auxI2 = zeros(size(pivData1.imPairNo));
        for ki = 1:numel(auxI1)
            aux = find(pivData2.imPairNo == pivData1.imPairNo(ki),1);
            if numel(aux)>0
                auxI2(ki) = aux;
            else
                auxI1(ki) = 0;
            end
        end
        aux = ~logical(auxI1+auxI2==0);
        auxI1 = auxI1(aux);
        auxI2 = auxI2(aux);
        rebuildRequired = false;
        if size(pivData1.U,3) ~= size(pivData2.U,3), rebuildRequired = true; end
        if sum(logical(auxI1) ~= (1:numel(auxI1))') > 0, rebuildRequired = true; end
        if sum(logical(auxI2) ~= (1:numel(auxI2))') > 0, rebuildRequired = true; end
        if rebuildRequired
            aux1 = pivManipulateData('readTimeSlice',pivData1,1);
            aux2 = pivManipulateData('readTimeSlice',pivData2,1);
            aux1 = pivManipulateData('initSequenceData',aux1,numel(auxI1));
            aux2 = pivManipulateData('initSequenceData',aux2,numel(auxI2));
            for ki = 1:numel(auxI1)
                aux = pivManipulateData('readTimeSlice',pivData1,auxI1(ki));
                aux1 = pivManipulateData('writeTimeSlice',aux1,aux,ki);
                aux = pivManipulateData('readTimeSlice',pivData2,auxI2(ki));
                aux2 = pivManipulateData('writeTimeSlice',aux2,aux,ki);
            end
            pivData1 = aux1;
            pivData2 = aux2;
        end
        % pivData1 and pivData2 contains now time slices for the same times.
        % if multiplicator field is not present, add it
        if ~isfield(pivData1,'multiplicator')
            pivData1.multiplicator = ones(size(pivData1.U))*mult1;
        end
        if ~isfield(pivData2,'multiplicator')
            pivData2.multiplicator = ones(size(pivData2.U))*mult2;
        end
        % consider pivData1 as primary, and replace them by pivData2 if they have better score
        pivDataOut = pivData1;
        % check ccPeak fields for NaNs. If NaNs are present, replace them by almost zeros:
        auxNaN = isnan(pivData2.ccPeak);
        pivData2.ccPeak(auxNaN) = 0.001;
        auxNaN = isnan(pivDataOut.ccPeak);
        pivDataOut.ccPeak(auxNaN) = 0.001;
        % check ccPeakSecondary fields for NaNs. If NaNs are present, replace them by ccPeak:
        auxNaN = isnan(pivData2.ccPeakSecondary);
        pivData2.ccPeakSecondary(auxNaN) = pivData2.ccPeak(auxNaN);
        auxNaN = isnan(pivDataOut.ccPeakSecondary);
        pivDataOut.ccPeakSecondary(auxNaN) = pivDataOut.ccPeak(auxNaN);
        % check status of pivData2. If replaced, decrease corresponding score.
        auxReplaced = logical(bitget(pivData2.Status,5));
        pivData2.multiplicator(auxReplaced) = 0;
        % compute velocity magnitude and clip it
        Umag2 = sqrt(pivData2.U.^2 + pivData2.V.^2);
        Umag2(logical(Umag2>UmagMax)) = UmagMax;
        UmagOut = sqrt(pivDataOut.U.^2 + pivDataOut.V.^2);
        UmagOut(logical(UmagOut>UmagMax)) = UmagMax;
        % check, where data should be taken from pivData2
        scoreOut = max((pivDataOut.ccPeak-pivDataOut.ccPeakSecondary-minDetecteability),0);
        scoreOut = pivDataOut.multiplicator.* UmagOut.*scoreOut;
        score2 = max((pivData2.ccPeak-pivData2.ccPeakSecondary-minDetecteability),0);
        score2 =  pivData2.multiplicator.* Umag2.*score2;
        auxShorterIsBetter = logical(scoreOut < score2);
        pivDataOut.U(auxShorterIsBetter) = pivData2.U(auxShorterIsBetter);
        pivDataOut.V(auxShorterIsBetter) = pivData2.V(auxShorterIsBetter);
        pivDataOut.Status(auxShorterIsBetter) = pivData2.Status(auxShorterIsBetter);
        pivDataOut.ccPeak(auxShorterIsBetter) = pivData2.ccPeak(auxShorterIsBetter);
        pivDataOut.ccPeakSecondary(auxShorterIsBetter) = pivData2.ccPeakSecondary(auxShorterIsBetter);
        pivDataOut.multiplicator(auxShorterIsBetter) = pivData2.multiplicator(auxShorterIsBetter);
        
        
    % ACTION toSingle
    case 'tosingle'
        pivDataIn = varargin{1};
        pivDataOut = pivDataIn;
        pivDataOut.X = single(pivDataOut.X);
        pivDataOut.Y = single(pivDataOut.Y);
        pivDataOut.U = single(pivDataOut.U);
        pivDataOut.V = single(pivDataOut.V);
        pivDataOut.ccPeak = single(pivDataOut.ccPeak);
        pivDataOut.ccPeakSecondary = single(pivDataOut.ccPeakSecondary);
        pivDataOut.Status = uint16(pivDataOut.Status);
        
        
    % ACTION toDouble
    case 'todouble'
        pivDataIn = varargin{1};
        pivDataOut = pivDataIn;
        pivDataOut.X = double(pivDataOut.X);
        pivDataOut.Y = double(pivDataOut.Y);
        pivDataOut.U = double(pivDataOut.U);
        pivDataOut.V = double(pivDataOut.V);
        pivDataOut.ccPeak = double(pivDataOut.ccPeak);
        pivDataOut.ccPeakSecondary = double(pivDataOut.ccPeakSecondary);
        pivDataOut.Status = uint16(pivDataOut.Status);
        
        
        % ACTIONS limitX, limitY
    case {'limitx','limity','limitxy'}
        pivDataIn = varargin{1};
        % get limits
        switch lower(action)
            case 'limitx'
                xmin = varargin{2}(1);
                xmax = varargin{2}(2);
                ymin = -Inf;
                ymax = +Inf;
            case 'limity'
                ymin = varargin{2}(1);
                ymax = varargin{2}(2);
                xmin = -Inf;
                xmax = +Inf;
            case 'limitxy'
                xmin = varargin{2}(1);
                xmax = varargin{2}(2);
                ymin = varargin{2}(3);
                ymax = varargin{2}(4);
        end
        % get range for x
        if isfield(pivDataIn,'X')
            xkeep = logical(pivDataIn.X(1,:) >= xmin) & logical(pivDataIn.X(1,:) <= xmax);
            ykeep = logical(pivDataIn.Y(:,1) >= ymin) & logical(pivDataIn.Y(:,1) <= ymax);
        end
        if isfield(pivDataIn,'spX')
            spXkeep = logical(pivDataIn.spX(1,:) >= xmin) & logical(pivDataIn.spX(1,:) <= xmax);
            spYkeep = logical(pivDataIn.spY(:,1) >= ymin) & logical(pivDataIn.spY(:,1) <= ymax);
        end
        if isfield(pivDataIn,'spuriousX')
            spurKeep = logical(pivDataIn.spuriousX(:,1) >= xmin) & logical(pivDataIn.spuriousX(:,1) <= xmax) & ...
                logical(pivDataIn.spuriousY(:,1) >= ymin) & logical(pivDataIn.spuriousY(:,1) <= ymax);
        end
        if isfield(pivDataIn,'spSpuriousX')
            spSpurKeep = logical(pivDataIn.spSpuriousX(:,1) >= xmin) & logical(pivDataIn.spSpuriousX(:,1) <= xmax) & ...
                logical(pivDataIn.spSpuriousY(:,1) >= ymin) & logical(pivDataIn.spSpuriousY(:,1) <= ymax);
        end
        % restrict the data
        if isfield(pivDataIn,'X'),pivDataOut.X = pivDataIn.X(ykeep,xkeep); end;
        if isfield(pivDataIn,'Y'),pivDataOut.Y = pivDataIn.Y(ykeep,xkeep); end;
        if isfield(pivDataIn,'iaStepX'),pivDataOut.iaStepX = pivDataIn.iaStepX; end;
        if isfield(pivDataIn,'iaStepY'),pivDataOut.iaStepY = pivDataIn.iaStepY; end;
        if isfield(pivDataIn,'Umean'),pivDataOut.Umean = pivDataIn.Umean(ykeep,xkeep); end;
        if isfield(pivDataIn,'Vmean'),pivDataOut.Vmean = pivDataIn.Vmean(ykeep,xkeep); end;
        if isfield(pivDataIn,'k'),pivDataOut.k = pivDataIn.k(ykeep,xkeep); end;
        if isfield(pivDataIn,'N'),pivDataOut.N = numel(pivDataOut.X); end;
        if isfield(pivDataIn,'U') && numel(size(pivDataIn.U))==2,pivDataOut.U = pivDataIn.U(ykeep,xkeep); end;
        if isfield(pivDataIn,'V') && numel(size(pivDataIn.V))==2,pivDataOut.V = pivDataIn.V(ykeep,xkeep); end;
        if isfield(pivDataIn,'vorticity') && numel(size(pivDataIn.vorticity))==2,pivDataOut.vorticity = pivDataIn.vorticity(ykeep,xkeep); end;
        if isfield(pivDataIn,'epsLEPIV') && numel(size(pivDataIn.epsLEPIV))==2,pivDataOut.epsLEPIV = pivDataIn.epsLEPIV(ykeep,xkeep); end;
        if isfield(pivDataIn,'epsMeanLEPIV') && numel(size(pivDataIn.epsMeanLEPIV))==2,pivDataOut.epsMeanLEPIV = pivDataIn.epsMeanLEPIV(ykeep,xkeep); end;
        if isfield(pivDataIn,'Status') && numel(size(pivDataIn.Status))==2,pivDataOut.Status = pivDataIn.Status(ykeep,xkeep); end;        
        if isfield(pivDataIn,'ccPeak') && numel(size(pivDataIn.ccPeak))==2,pivDataOut.ccPeak = pivDataIn.ccPeak(ykeep,xkeep); end;        
        if isfield(pivDataIn,'ccPeakSecondary') && numel(size(pivDataIn.ccPeakSecondary))==2,pivDataOut.ccPeakSecondary = pivDataIn.ccPeakSecondary(ykeep,xkeep); end;        
        if isfield(pivDataIn,'U') && numel(size(pivDataIn.U))==3,pivDataOut.U = pivDataIn.U(ykeep,xkeep,:); end;
        if isfield(pivDataIn,'V') && numel(size(pivDataIn.V))==3,pivDataOut.V = pivDataIn.V(ykeep,xkeep,:); end;
        if isfield(pivDataIn,'vorticity') && numel(size(pivDataIn.vorticity))==3,pivDataOut.vorticity = pivDataIn.vorticity(ykeep,xkeep,:); end;
        if isfield(pivDataIn,'epsLEPIV') && numel(size(pivDataIn.epsLEPIV))==3,pivDataOut.epsLEPIV = pivDataIn.epsLEPIV(ykeep,xkeep,:); end;
        if isfield(pivDataIn,'Status') && numel(size(pivDataIn.Status))==3,pivDataOut.Status = pivDataIn.Status(ykeep,xkeep,:); end;
        if isfield(pivDataIn,'ccPeak') && numel(size(pivDataIn.ccPeak))==3,pivDataOut.ccPeak = pivDataIn.ccPeak(ykeep,xkeep,:); end;
        if isfield(pivDataIn,'ccPeakSecondary') && numel(size(pivDataIn.ccPeakSecondary))==3,pivDataOut.ccPeakSecondary = pivDataIn.ccPeakSecondary(ykeep,xkeep,:); end;
        if isfield(pivDataIn,'spuriousX'),pivDataOut.spuriousX = pivDataIn.spuriousX(spurKeep); end;
        if isfield(pivDataIn,'spuriousY'),pivDataOut.spuriousY = pivDataIn.spuriousY(spurKeep); end;
        if isfield(pivDataIn,'spuriousU'),pivDataOut.spuriousU = pivDataIn.spuriousU(spurKeep); end;
        if isfield(pivDataIn,'spuriousV'),pivDataOut.spuriousV = pivDataIn.spuriousV(spurKeep); end;
        if isfield(pivDataIn,'spuriousX'),pivDataOut.spuriousN = numel(pivDataOut.spuriousX); end
        if isfield(pivDataIn,'imFilename1'),pivDataOut.imFilename1 = pivDataIn.imFilename1; end;
        if isfield(pivDataIn,'imFilename2'),pivDataOut.imFilename2 = pivDataIn.imFilename2; end;
        if isfield(pivDataIn,'imMaskFilename1'),pivDataOut.imMaskFilename1 = pivDataIn.imMaskFilename1; end;
        if isfield(pivDataIn,'imMaskFilename2'),pivDataOut.imMaskFilename2 = pivDataIn.imMaskFilename2; end;
        if isfield(pivDataIn,'imNo1'),pivDataOut.imNo1 = pivDataIn.imNo1; end;
        if isfield(pivDataIn,'imNo2'),pivDataOut.imNo2 = pivDataIn.imNo2; end;
        if isfield(pivDataIn,'spCC'),pivDataOut.spCC = pivDataIn.spCC(spYkeep,spXkeep,:,:); end;        
        if isfield(pivDataIn,'spAC'),pivDataOut.spAC = pivDataIn.spAC(spYkeep,spXkeep,:,:); end;   
        if isfield(pivDataIn,'spX'),pivDataOut.spX = pivDataIn.spX(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spY'),pivDataOut.spY = pivDataIn.spY(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spStatus'),pivDataOut.spStatus = pivDataIn.spStatus(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spU0'),pivDataOut.spU0 = pivDataIn.spU0(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spV0'),pivDataOut.spV0 = pivDataIn.spV0(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spCCmax'),pivDataOut.spCCmax = pivDataIn.spCCmax(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spUint'),pivDataOut.spUint = pivDataIn.spUint(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spVint'),pivDataOut.spVint = pivDataIn.spVint(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spC0int'),pivDataOut.spC0int = pivDataIn.spC0int(spYkeep,spXkeep); end;  
        if isfield(pivDataIn,'spC1int'),pivDataOut.spC1int = pivDataIn.spC1int(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spC2int'),pivDataOut.spC2int = pivDataIn.spC2int(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spPhiint'),pivDataOut.spPhiint = pivDataIn.spPhiint(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spUfit'),pivDataOut.spUfit = pivDataIn.spUfit(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spVfit'),pivDataOut.spVfit = pivDataIn.spVfit(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spC0fit'),pivDataOut.spC0fit = pivDataIn.spC0fit(spYkeep,spXkeep); end;          
        if isfield(pivDataIn,'spC1fit'),pivDataOut.spC1fit = pivDataIn.spC1fit(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spC2fit'),pivDataOut.spC2fit = pivDataIn.spC2fit(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spPhifit'),pivDataOut.spPhifit = pivDataIn.spPhifit(spYkeep,spXkeep); end;   
        if isfield(pivDataIn,'spSpuriousX'),pivDataOut.spSpuriousX = pivDataIn.spSpuriousX(spSpurKeep); end;   
        if isfield(pivDataIn,'spSpuriousY'),pivDataOut.spSpuriousY = pivDataIn.spSpuriousY(spSpurKeep); end;   
        if isfield(pivDataIn,'spSpuriousU'),pivDataOut.spSpuriousU = pivDataIn.spSpuriousU(spSpurKeep); end;   
        if isfield(pivDataIn,'spSpuriousV'),pivDataOut.spSpuriousV = pivDataIn.spSpuriousV(spSpurKeep); end;   
        if isfield(pivDataIn,'spSpuriousCx'),pivDataOut.spSpuriousCx = pivDataIn.spSpuriousCx(spSpurKeep); end;
        if isfield(pivDataIn,'spSpuriousCy'),pivDataOut.spSpuriousCy = pivDataIn.spSpuriousCy(spSpurKeep); end;
        if isfield(pivDataIn,'spSpuriousPhi'),pivDataOut.spSpuriousPhi = pivDataIn.spSpuriousPhi(spSpurKeep); end;
        if isfield(pivDataIn,'spSpuriousCcmax'),pivDataOut.spSpuriousCcmax = pivDataIn.spSpuriousCcmax(spSpurKeep); end;
        if isfield(pivDataIn,'spDeltaXNeg'),pivDataOut.spDeltaXNeg = pivDataIn.spDeltaXNeg; end;   
        if isfield(pivDataIn,'spDeltaXPos'),pivDataOut.spDeltaXPos = pivDataIn.spDeltaXPos; end;   
        if isfield(pivDataIn,'spDeltaYNeg'),pivDataOut.spDeltaYNeg = pivDataIn.spDeltaYNeg; end;   
        if isfield(pivDataIn,'spDeltaYPos'),pivDataOut.spDeltaYPos = pivDataIn.spDeltaYPos; end;   
        if isfield(pivDataIn,'spBindX'),pivDataOut.spBindX = pivDataIn.spBindX; end;   
        if isfield(pivDataIn,'spBindY'),pivDataOut.spBindY = pivDataIn.spBindY; end;   
        % copy all fields starting with 'exp'
        fields = fieldnames(pivDataIn);
        auxI = strfind(fields,'exp');
        for kk = 1:numel(auxI)
            if numel(auxI{kk})>0 && auxI{kk,1}==1
                pivDataOut.(fields{kk}) = pivDataIn.(fields{kk});
            end
        end
        % copy all fields starting with 'im'
        fields = fieldnames(pivDataIn);
        auxI = strfind(fields,'im');
        for kk = 1:numel(auxI)
            if numel(auxI{kk})>0 && auxI{kk,1}==1
                pivDataOut.(fields{kk}) = pivDataIn.(fields{kk});
            end
        end
        
    % ACTION not recognized
    otherwise
        fprintf('Error: action %s is not recognized by pivManipulateData.\n',action);
end
end


%% LOCAL FUNCTIONS

function [D] = diffField(F,X,dim)
% computes central difference dF/dX, where F is 2D or 3D field, X is a 2D
% matrix with coordinates; dim defines, along which dimension is X varied
% (use dim = 1 for differentiating across rows (i.e. dF/dY in PIVsuite),
% and dim=2 for differentiating across columns (dF/dX)).

% initialize fields
F1 = zeros(size(F,1)+2,size(F,2)+2,size(F,3)) + NaN;
F2 = F1;
X1 = zeros(size(X,1)+2,size(X,2)+2) + NaN;
X2 = X1;
% create shifted and padded arrays
if dim==1   % derivative along Y
    F2(3:end,  2:end-1,:) = F;
    F1(1:end-2,2:end-1,:) = F;
    X2(3:end,  2:end-1,:) = X;
    X1(1:end-2,2:end-1,:) = X;
elseif dim==2
    F2(2:end-1,3:end,  :) = F;
    F1(2:end-1,1:end-2,:) = F;
    X2(2:end-1,3:end,  :) = X;
    X1(2:end-1,1:end-2,:) = X;
else
    error('pivManipulateData.m:diffField: invalid value of dim');
end
% differentiate
D = F1 + NaN;
for kk=1:size(F,3)
    D(:,:,kk) = (F2(:,:,kk)-F1(:,:,kk))./(X2(:,:)-X1(:,:));
end
% clip padded values
D = D(2:end-1,2:end-1,:);
end