piv_3d_10.m

function piv_3d_10(piv_parameters);
% This function is designed to process two or three dimensional PIV data using
% the parameters specified by the 'piv_parameters' data structure.
%
% This code is based upon the code 'basic_3d_rpc_processing_05.m'.
%
% Updates on previous versions:
%
%   Version 08
%
%       Saves the PIV processing parameters data structure to the PIV
%       vector field write directory for later reference.
%
%   Version 09
%
%       Added the ability to use several different outlier vector
%       replacement methods including a local mean interpolation (default
%       in previous versions), a Laplacian interpolation, and a Delaunay
%       interpolation.
%
%   Version 10
%
%       This updated the order of the spatial and temporal dimensions to be
%       consistent with other PIV software versions and to minimize the
%       number of overhead operations (since MATLAB prefers some orders
%       over others).
%
%       The ability to import other data formats was also included.  The
%       code can now import 2D images in one of the following formats: 
%       'bmp', 'jpg', 'jp2', 'png', 'ppm', 'tif', 'mat', 'dcm', 'cdf', or 
%       'h5' while 3D images may be one of the following formats: 'mat', 
%       'dcm', 'cdf', or 'h5'.
%
%       Any NaN velocity vector values that are produced during the 
%       correlations (likely due to simulated data not being seeded within
%       a window) are considered outliers and replaced accordingly.
%
% Things to work on, modify, fix, et cetera:
%
%   Change the code that is optimized for C to code that is optimized for
%   matlab.
%
%   Check on the speed of the if-then statements added to make the code run
%   for 2D images as well.
%
%   Add peak fitting algorithms besides the 3-point Gaussian.
%
%   Add the ability to record multiple correlation peaks, specifically to
%   be used during validation.  Possibly also record the peak ratio.
%
%   Add the ability to do multiple iterations of validation and replacement
%   for each pass.
%
%   Add in the ability to perform deform correlations.
%
% Authors:		Rod La Foy
% First Written On:	31 October 2014
% Last Modified On:	21 November 2014

% This extracts the vector field saving directory from the data processing structure
vector_field_write_directory=piv_parameters.general_parameters.vector_field_write_directory;
% This saves a copy of the processing parameters to the vector field write
% directory for future reference
save([vector_field_write_directory,'piv_processing_parameters.mat'],'piv_parameters');

% This extracts the number of passes to perform
pass_number=piv_parameters.general_parameters.pass_number;

% This is the initial image frame to load
initial_image_frame=piv_parameters.general_parameters.initial_image_frame;
% This is the final image frame to load
final_image_frame=piv_parameters.general_parameters.final_image_frame;
% This is the number of frames to step each iteration
image_frame_step=piv_parameters.general_parameters.image_frame_step;
% This is the number of frames to step between correlation frame pairs
image_correlation_step=piv_parameters.general_parameters.image_correlation_step;

% This is the Boolean value stating whether to perform window deformation
perform_window_deformation=piv_parameters.general_parameters.perform_window_deformation;
% If window deformation is to be performed, this loads the window
% deformation parameters
if perform_window_deformation;
    % This is the minimum number of window deformation operations to
    % perform
    window_deformation_iteration_min=piv_parameters.general_parameters.window_deformation_iteration_min;
    % This is the maximum number of window deformation operations to
    % perform
    window_deformation_iteration_max=piv_parameters.general_parameters.window_deformation_iteration_max;
    % This is the convergence threshhold of the window deformation process
    window_deformation_threshhold=piv_parameters.general_parameters.window_deformation_threshhold;
end;

% This is a Boolean value stating whether to perform pyramid correlations
perform_pyramid_correlations=piv_parameters.general_parameters.perform_pyramid_correlations;

% This is the image filename extension
image_filename_extension=piv_parameters.general_parameters.image_filename_extension;
% This is the image variable name (if the variable is saved within a data
% object that contain multiple variables, ie a 'mat' file)
image_variable_name=piv_parameters.general_parameters.image_variable_name;

% This is a list of the images within the image read directory
image_filename_list=dir([piv_parameters.general_parameters.image_read_directory,piv_parameters.general_parameters.image_filename_prefix,'*',image_filename_extension]);

% This is the filename of the first image file (for the purpose of determing the image
% file size)
first_image_filename=[piv_parameters.general_parameters.image_read_directory,image_filename_list(1).name];

% This finds the size of the images being loaded which can be 'bmp', 'jpg',
% 'jp2', 'png', 'ppm', or 'tif' images in 2D or 'mat' or 'dcm' images in 3D
if any(strcmpi(image_filename_extension,{'bmp','jpg','jp2','png','ppm','tif'}));
    
    % This loads in the first image
    I_First_Image=imread(first_image_filename);
    % This measures the image size
    image_size=size(I_First_Image);
    % This clears the image from memory
    clear('I_First_Image');
    % This adds the third dimension as 1 if it does not exist
    if length(image_size)==2;
        % This adds the third dimension of the image size as 1
        image_size(3)=1;
    end;
    % This clears the first image from memory
    clear('I_First_Image');
    
elseif strcmpi(image_filename_extension,'mat');
    
    % This creates a matlab mat object containing information about the image
    first_image_object=matfile(first_image_filename);
    % This extracts the image size information
    image_size=size(first_image_object,image_variable_name);
    % This adds the third dimension as 1 if it does not exist
    if length(image_size)==2;
        % This adds the third dimension of the image size as 1
        image_size(3)=1;
    end;
    
elseif strcmpi(image_filename_extension,'dcm');
    
    % This reads int the DICOM image
    I_First_Image=dicomread(first_image_filename);
    % This extracts the image size information
    image_size=size(I_First_Image);
    % This adds the third dimension as 1 if it does not exist and removes
    % the third dimension if the image is 3D since the 3rd dimension will
    % correspond to the colormap
    if length(image_size)==2;
        % This adds the third dimension of the image size as 1
        image_size(3)=1;
    elseif length(image_size)==4;
        % This removes the 3rd dimension since this corresponds to the
        % colormap fod DICOM files which is not used for PIV data
        image_size(3)=[];
    end;
    % This clears the first image from memory
    clear('I_First_Image');
    
elseif strcmpi(image_filename_extension,'cdf');
    
    % This reads in the Common Data Format image into a cell array
    I_First_Image=cdfread(first_image_filename,'Variables',image_variable_name);
    % This extracts the image size information from the cell array
    image_size=size(I_First_Image{1});
    % This adds the third dimension as 1 if it does not exist and removes
    % the third dimension if the image is 3D since the 3rd dimension will
    % correspond to the colormap
    if length(image_size)==2;
        % This adds the third dimension of the image size as 1
        image_size(3)=1;
    end;
    % This clears the first image from memory
    clear('I_First_Image');
    
elseif strcmpi(image_filename_extension,'h5');
    
    % This reads in the HDF5 image
    I_First_Image=hdf5read(first_image_filename,image_variable_name);
    % This extracts the image size information
    image_size=size(I_First_Image);
    % This adds the third dimension as 1 if it does not exist and removes
    % the third dimension if the image is 3D since the 3rd dimension will
    % correspond to the colormap
    if length(image_size)==2;
        % This adds the third dimension of the image size as 1
        image_size(3)=1;
    end;
    % This clears the first image from memory
    clear('I_First_Image');
    
end;

% These are the coordinates of the estimated velocities (initially nothing is assummed
% about the velocity field so a zero displacement is used and since the extrapolation value
% is set to 0, only one arbitrary coordinate must be specified)
if image_size(3)==1;
    % This calculates the initial interpolation coordinates with only one
    % kk coordinate value
    [jj_position,ii_position,kk_position]=meshgrid([1,image_size(1)],[1,image_size(2)],[1]);
    % This is the estimated (initially zero) velocity for the 2D array
    ii_velocity=zeros(2,2,1);
    jj_velocity=zeros(2,2,1);
    kk_velocity=zeros(2,2,1);
    
    
elseif image_size(3)>1;
    % This calculates the initial interpolation coordinates with the full
    % 3D array of coordinates
    [jj_position,ii_position,kk_position]=meshgrid([1,image_size(1)],[1,image_size(2)],[1,image_size(3)]);
    % This is the estimated (initially zero) velocity for the 3D array
    ii_velocity=zeros(1,1,1);   % originally 2,2,2
    jj_velocity=zeros(1,1,1);   % originally 2,2,2
    kk_velocity=zeros(1,1,1);   % originally 2,2,2
    
    
end;

% This iterates through the passes
for pass_index=1:pass_number;
    
    % This displays that the current pass is being processed
    fprintf('\n\nCompleting pass %d of %d passes.\n',pass_index,pass_number);
    
    % This extracts the effective current window resolution
    window_resolution=piv_parameters.pass_parameters(pass_index).window_resolution;
    % This extracts the full window size
    window_size=piv_parameters.pass_parameters(pass_index).window_size;
    % This extracts the window gridding method for the current pass
    window_gridding_method=piv_parameters.pass_parameters(pass_index).window_gridding_method;
    
    % This is the Boolean value stating whether to perform validation on the vector field
    validate_vector_field=piv_parameters.pass_parameters(pass_index).validate_vector_field;
    
    % This is the Boolean value stating whether to perform smoothing of the velocity vector field
    smooth_vector_field=piv_parameters.pass_parameters(pass_index).smooth_vector_field;
    
    % If specified by the user to perform smoothing of the velocity field, this loads the
    % smoothing specific parameters
    if smooth_vector_field;
        
        % This loads the Gaussian smoothing standard deviation
        gaussian_smoothing_kernal_std=piv_parameters.pass_parameters(pass_index).gaussian_smoothing_kernal_std;
        % This loads the Gaussian smoothing kernal size
        gaussian_smoothing_kernal_size=piv_parameters.pass_parameters(pass_index).gaussian_smoothing_kernal_size;
        
    end;
    
    % If specified by the user to perform validation, this loads the loads the validation
    % specific parameters
    if validate_vector_field;
        
        % This loads the vector outlier threshhold limits
        minimum_outlier_threshhold=piv_parameters.pass_parameters(pass_index).minimum_outlier_threshhold;
        maximum_outlier_threshhold=piv_parameters.pass_parameters(pass_index).maximum_outlier_threshhold;
        
        % This loads the UOD kernal size
        uod_kernal_size=piv_parameters.pass_parameters(pass_index).uod_kernal_size;
        % This loads the UOD expected error
        uod_epsilon=piv_parameters.pass_parameters(pass_index).uod_epsilon;
        % This loads the UOD residual threshhold value
        uod_residual_threshhold=piv_parameters.pass_parameters(pass_index).uod_residual_threshhold;
        
        % This loads the outlier vector replacement method
        vector_replacement_method=piv_parameters.pass_parameters(pass_index).vector_replacement_method;
        % This loads the relevant vector replacement parameters based upon
        % the interpolation method
        if strcmp(vector_replacement_method,'local_mean');
            % This loads the minimum number of local vectors used to
            % calculate the local mean
            local_mean_minimum_valid_vector_number=piv_parameters.pass_parameters(pass_index).local_mean_minimum_valid_vector_number;
        elseif strcmp(vector_replacement_method,'laplacian_interpolation');
            % This loads the number of adjacent points to use in 
            % calculating the interpolated values (this can be either 4 or 
            % 8 in 2D and 6, 18, or 26 in 3D)
            laplacian_interpolation_adjacent_connectivity=piv_parameters.pass_parameters(pass_index).laplacian_interpolation_adjacent_connectivity;
        elseif strcmp(vector_replacement_method,'delaunay_interpolation');
            % This loads the Delaunay interpolation weighting method to be
            % used in interpolating values
            delaunay_interpolation_weighting_method=piv_parameters.pass_parameters(pass_index).delaunay_interpolation_weighting_method;
        end;
        
    end;
    
    % This calculates the grid spacing for the current pass using either
    % the amount of window overlap or the window spacing
    if strcmp(window_gridding_method,'window_overlap');
        % This is the ratio of the window overlap distance
        window_overlap=piv_parameters.pass_parameters(pass_index).window_overlap;
        % This is the spacing between the centers of the windows
        grid_spacing=round(window_resolution.*(1-window_overlap));  %originally grid_spacing=round(window_resolution.*(1-window_overlap));
        % This checks if the overlap percentage is so high that grid spacing is zero (although
        % this should actually never be done)
        if grid_spacing==0;
            % This sets the grid spacing to 1 voxel
            grid_spacing=ones(size(grid_spacing));
        end;
    elseif strcmp(window_gridding_method,'window_spacing');
        % This extracts the window grid spacing from the parameters structure
        grid_spacing=piv_parameters.pass_parameters(pass_index).window_spacing;
    end;
    

    
    
    % This calculates the window domains
    [window_min,window_max,window_center,window_number]=calculate_window_domains(image_size,window_resolution,window_size,grid_spacing);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This iterates through the frames calculating the velocity fields.                      %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Z_ensemble = zeros(window_size);
    
    % This iterates through the image frames calculating the vector fields
    for frame_index=initial_image_frame:image_frame_step:final_image_frame;
        
        % This displays that the current frame is being processed
        fprintf('\nProcessing frame %d of %d frames.\n',frame_index,final_image_frame);
        
        
        
        % If this is the second or higher pass, this loads the
        % previously calculated vector field which is used to deform
        % the images in the pyramid correlations
        if pass_index>1;
            
            % This is the filename to load the data from
            data_filename_read=[vector_field_write_directory,'pass_',sprintf('%02.0f',pass_index-1),'_frame_',sprintf('%04.0f',frame_index),'_x_',sprintf('%04.0f',frame_index+image_correlation_step),'.mat'];
            % This loads in the data file
            load(data_filename_read);
            
            % This renames the data loaded from memory to be conistent 
            % with the variables names within this code (ie Prana uses X,
            % Y, Z, U, V, W for the position and velocity values and this
            % code uses the index dimensions ii, jj, kk)
            ii_position=Y;
            jj_position=X;
            kk_position=Z;
            ii_velocity=V;
            jj_velocity=U;
            kk_velocity=W;
            
        end;
        

        
        % This checks whether standard correlations are supposed to be
        % performed
        if not(perform_pyramid_correlations);
            
            % This performs the standard correlations of the current frame
            [Z,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=standard_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity);
            
            %ensemble correlation addition
            Z_ensemble = Z_ensemble + Z;
            C=ifftn(Z_ensemble,'symmetric');
            
            displacement_vector=subpixel_displacement_vector(C,window_size);
            
            
            
        elseif perform_pyramid_correlations;
            
            % This performs the pyramid correlations of the current frame
            [ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=pyramid_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity);
            
        end;
        
        
        
        
        
        % This performs validation of the vector field if specified by the user
        if validate_vector_field;
            
            % This displays that the vector field is being validated
            disp('Validating the vector field . . . ');
            
            % This determines the outliers of the vector field based upon several metrics
            outlier_vector_array=locate_vector_outliers(ii_velocity,jj_velocity,kk_velocity,minimum_outlier_threshhold,maximum_outlier_threshhold,uod_kernal_size,uod_epsilon,uod_residual_threshhold);
            
            % This replaces the outlier vectors using the specified method
            if strcmp(vector_replacement_method,'local_mean');
                % This replaces the outlier vectors using a local mean of valid vectors
                [ii_velocity,jj_velocity,kk_velocity]=local_mean_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,local_mean_minimum_valid_vector_number);
            elseif strcmp(vector_replacement_method,'laplacian_interpolation');
                % This replaces the outlier vectors using Laplacian interpolation
                [ii_velocity,jj_velocity,kk_velocity]=laplacian_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,laplacian_interpolation_adjacent_connectivity);
            elseif strcmp(vector_replacement_method,'delaunay_interpolation');
                % This replaces the outlier vectors using a weighted sum based upon the
                % Delaunay triangulation of the valid vectors
                [ii_velocity,jj_velocity,kk_velocity]=delaunay_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,delaunay_interpolation_weighting_method);
            end;
            
        end;
        
        % This performs smoothing of the vector field if specified by the users
        if smooth_vector_field;
            
            % This displays that the vector field is being smoothed
            disp('Smoothing the vector field . . . ');
            
            % This smooths the vector fields
            [ii_velocity,jj_velocity,kk_velocity]=smooth_velocity_field(ii_velocity,jj_velocity,kk_velocity,gaussian_smoothing_kernal_size,gaussian_smoothing_kernal_std*ones(1,3));
            
        end;
        
        % This is the filename to save the data as
        data_filename_write=[vector_field_write_directory,'pass_',sprintf('%02.0f',pass_index),'_frame_',sprintf('%04.0f',frame_index),'_x_',sprintf('%04.0f',frame_index+image_correlation_step),'.mat'];

        % This renames the data from the current vector field to be 
        % conistent with the variables names used by other vector field 
        % code (ie Prana uses X, Y, Z, U, V, W for the position and 
        % velocity values and this code uses the index dimensions ii, jj, 
        % kk)
        Y=ii_position;
        X=jj_position;
        Z=kk_position;
        V=ii_velocity;
        U=jj_velocity;
        W=kk_velocity;
        
        % This renames the validation array to be consistent with the
        % variables used by other vector field code (ie Prana uses Valid as
        % the name for the valid vectors and this code uses
        % outlier_vector_array)
        if validate_vector_field;
            % This creates the binary array variable 'Valid' which gives
            % the locations of the vectors that were considered valid
            Valid=not(outlier_vector_array);
        else;
            % This creates an array of NaNs for the variable 'Valid' since
            % validation was not performed on the vector field
            Valid=NaN(size(X));
        end;

        
        % This saves the vector field data to memory
        save(data_filename_write,'X','Y','Z','U','V','W','Valid','-v7.3');
        
    end;
    
end;



function [Z,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=standard_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity);
% This function calculates standard pair correlations of the data set listed
% within the parameters variable 'piv_parameters'.  The standard correlation
% consists of taking pair-wise correlations from a time series of images.
% The input parameters of this function are given by
%
%   piv_parameters      A data structure containing all the necessary
%                       parameters required to perform the PIV processing
%
%   pass_index          A scalar integer giving the current pass of the PIV
%                       processing
%
%   frame_index         A scalar integer giving the index of the current
%                       frame being processed from the directory of images.
%                       For pyramid correlations, the velocity field will
%                       always be calculated for the frame_index and
%                       frame_index + 1 froms, ie the changing the
%                       correlation step does nothing.
%
%   window_number       A 1 x 3 vector giving the number of PIV windows to
%                       process in the first, second, and third dimensions
%                       respectively.  If the images are 2D then
%                       window_number(3) should always equal 1.
%
%   window_min          A L x 3 vector giving the minimum indices of the
%                       window domains in each respective dimension where L
%                       is the total number of windows to process.
%
%   window_max          A L x 3 vector giving the maximum indices of the
%                       window domains in each respective dimension where L
%                       is the total number of windows to process.
%
%   window_center       A L x 3 vector giving the coordinates of the center
%                       of the window domains in each respective dimension
%                       where L is the total number of windows to process.
%
%   ii_position         A M x N x P array giving the 1st dimension
%                       coordinates of the previously measured vector field
%                       if this is at least the second pass.  If this is
%                       the first pass then the coordinates correspond to
%                       an all-zero vector field.
%
%   jj_position         A M x N x P array giving the 2nd dimension
%                       coordinates of the previously measured vector field
%                       if this is at least the second pass.  If this is
%                       the first pass then the coordinates correspond to
%                       an all-zero vector field.
%
%   kk_position         A M x N x P array giving the 3rd dimension
%                       coordinates of the previously measured vector field
%                       if this is at least the second pass.  If this is
%                       the first pass then the coordinates correspond to
%                       an all-zero vector field.
%
%   ii_velocity         A M x N x P array giving the 1st dimension velocity
%                       of the previously measured vector field if this is
%                       at least the second pass.  If this is the first
%                       pass then the velocity will be set to all zeros.
%
%   jj_velocity         A M x N x P array giving the 2nd dimension velocity
%                       of the previously measured vector field if this is
%                       at least the second pass.  If this is the first
%                       pass then the velocity will be set to all zeros.
%
%   kk_velocity         A M x N x P array giving the 3rd dimension velocity
%                       of the previously measured vector field if this is
%                       at least the second pass.  If this is the first
%                       pass then the velocity will be set to all zeros.
%
% The output parameters give the measured (un-validated, un-smoothed)
% velocity field of the current frame pair and are described by
%
%   ii_position         A Q x R x S array giving the 1st dimension
%                       coordinates of the measured vector field, ie the
%                       center of the correlation window in the 1st
%                       dimension.
%
%   jj_position         A Q x R x S array giving the 2nd dimension
%                       coordinates of the measured vector field, ie the
%                       center of the correlation window in the 1st
%                       dimension.
%
%   kk_position         A Q x R x S array giving the 3rd dimension
%                       coordinates of the measured vector field, ie the
%                       center of the correlation window in the 1st
%                       dimension.
%
%   ii_velocity         A Q x R x S array giving the 1st dimension velocity
%                       of the current vector field measured at each
%                       correlation window.
%
%   jj_velocity         A Q x R x S array giving the 1st dimension velocity
%                       of the current vector field measured at each
%                       correlation window.
%
%   kk_velocity         A Q x R x S array giving the 1st dimension velocity
%                       of the current vector field measured at each
%                       correlation window.
%
% Authors: Rod La Foy
% First Created On: 14 March 2014
% Last Modified On: 14 March 2014

% This is the image filename extension
image_filename_extension=piv_parameters.general_parameters.image_filename_extension;
% This is the image variable name (if the variable is saved within a data
% object that contain multiple variables, ie a 'mat' file)
image_variable_name=piv_parameters.general_parameters.image_variable_name;

% This is a list of the images within the image read directory
image_filename_list=dir([piv_parameters.general_parameters.image_read_directory,piv_parameters.general_parameters.image_filename_prefix,'*',image_filename_extension]);

% This extracts the effective current window resolution
window_resolution=piv_parameters.pass_parameters(pass_index).window_resolution;
% This extracts the full window size
window_size=piv_parameters.pass_parameters(pass_index).window_size;

% This is the number of frames to step between correlation frame pairs
image_correlation_step=piv_parameters.general_parameters.image_correlation_step;

% This is the bulk window offset distance
bulk_window_offset=piv_parameters.pass_parameters(pass_index).bulk_window_offset;

% This is the correlation method to use for the current pass
correlation_method=piv_parameters.pass_parameters(pass_index).correlation_method;
% If the correlation method is the RPC method, this creates the RPC filter
if strcmp(correlation_method,'RPC');
    % This extracts the size of the particles in the images
    particle_diameter=piv_parameters.general_parameters.particle_diameter;
    % This creates the spectral energy filter
    spectral_filter=single(fftshift(create_spectral_filter(window_size,particle_diameter)));
elseif strcmp(correlation_method,'SCC');
    % This sets the spectral filter to a null value
    spectral_filter=[];
end;

% This is a Boolean value stating whether to zero-mean the correlation
% windows
zero_mean_windows=piv_parameters.pass_parameters(pass_index).zero_mean_windows;

% This sets the fft library for optimal speed calculation
fftw('planner','measure');

% This determines the size of the Gaussian function to use for the window masking
gaussian_width=determine_gaussian_size_2(window_size,window_resolution);

% This creates a Gaussian windowing function
gaussian_filter=single(create_gaussian_filter(window_size,gaussian_width));

% These are the coordinate vectors
ii_position_current=reshape(window_center(:,1),window_number);
jj_position_current=reshape(window_center(:,2),window_number);
kk_position_current=reshape(window_center(:,3),window_number);


% This initializes the velocity vectors
ii_velocity_current=zeros(size(ii_position_current));
jj_velocity_current=zeros(size(ii_position_current));
kk_velocity_current=zeros(size(ii_position_current));

% This determines the domain over which the known velocity field should be
% extrapolated
% ii_extrap_min=min(ii_position_current(:))-ceil(window_size(1)/2);
% ii_extrap_max=max(ii_position_current(:))+ceil(window_size(1)/2);
% jj_extrap_min=min(jj_position_current(:))-ceil(window_size(2)/2);
% jj_extrap_max=max(jj_position_current(:))+ceil(window_size(2)/2);
% kk_extrap_min=min(kk_position_current(:))-ceil(window_size(3)/2);
% kk_extrap_max=max(kk_position_current(:))+ceil(window_size(3)/2);
% % This adds extrapolated points to the previously measured vector field so
% % that if the vector field is interpolated at points near the edge of the
% % image, the interpolated points will have approximately correct values (as
% % opposed to zero or NaNs)
% [~,~,~,ii_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,ii_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max);
% [~,~,~,jj_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,jj_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max);
% [ii_position,jj_position,kk_position,kk_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,kk_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max);

% If the bulk window offset is non-zero, this adds it to the previous pass
% velocity field estimate
if any(bulk_window_offset);
    % This adds the bulk window offset to the previously measured velocity
    % field
    ii_velocity=ii_velocity+bulk_window_offset(1);
    jj_velocity=jj_velocity+bulk_window_offset(2);
    kk_velocity=kk_velocity+bulk_window_offset(3);
    % This creates a Boolean value stating to apply the bulk window offset
    perform_bulk_window_offset=true;
else;
    % This creates a Boolean value stating to apply the bulk window offset
    perform_bulk_window_offset=false;
end;

% If this is the second pass or higher, this calculates the discrete window
% offset, otherwise the offset is just set to zero
if (pass_index==1)&&(not(perform_bulk_window_offset));
    % This sets the first frame minimum and maximum values to the
    % non-offset values
    frame_1_window_min=window_min;
    frame_1_window_max=window_max;
    % This sets the second frame minimum and maximum values to the
    % non-offset values
    frame_2_window_min=window_min;
    frame_2_window_max=window_max;
    % This sets the displacement offset array to all zeros
    displacement_offset=zeros(size(window_min));
else;
    % This function adds the known velocity field to the window domains, ie
    % accounts for the discrete window offset
    [frame_1_window_min,frame_1_window_max,frame_2_window_min,frame_2_window_max,displacement_offset]=dwo_window_domains(window_min,window_max,window_center,image_correlation_step,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity);
end;

% This is the filename of the first image to load
image_1_filename=[piv_parameters.general_parameters.image_read_directory,image_filename_list(frame_index).name];
% This is the filename of the first image to load


image_2_filename=[piv_parameters.general_parameters.image_read_directory,image_filename_list(frame_index+image_correlation_step).name];

% This finds the size of the images being loaded which can be 'bmp', 'jpg',
% 'jp2', 'png', 'ppm', or 'tif' images in 2D or 'mat images in 3D
if any(strcmpi(image_filename_extension,{'bmp','jpg','jp2','png','ppm','tif'}));
    
    % This loads the first image
    I1=double(imread(image_1_filename));
    % This loads the second image
    I2=double(imread(image_2_filename));
    
elseif strcmpi(image_filename_extension,'mat');
    
    % This creates a matlab file object of the first image of the correlation pairs
    image_1_object=matfile(image_1_filename);
    % This loads the first image
    eval(['I1=double(image_1_object.',image_variable_name,');']);
    % This creates a matlab file object of the second image of the correlation pairs
    image_2_object=matfile(image_2_filename);
    % This loads the second image
    eval(['I2=double(image_2_object.',image_variable_name,');']);
    
elseif strcmpi(image_filename_extension,'dcm');
    
    % This loads the first image
    I1=double(squeeze(dicomread(image_1_filename)));
    % This loads the second image
    I2=double(squeeze(dicomread(image_2_filename)));
    
elseif strcmpi(image_filename_extension,'cdf');
    
    % This loads the first image cell array
    I1_Temp=double(cdfread(image_1_filename,'Variables',image_variable_name));
    % This loads the image from the cell array
    I1=I1_Temp{1};
    % This clears the cell array from memory
    clear('I1_Temp');
    % This loads the second image
    I2_Temp=double(cdfread(image_2_filename,'Variables',image_variable_name));
    % This loads the image from the cell array
    I2=I2_Temp{1};
    % This clears the cell array from memory
    clear('I2_Temp');
    
elseif strcmpi(image_filename_extension,'h5');
    
    % This loads the first image
    I1=double(hdf5read(image_1_filename,image_variable_name));
    % This loads the second image
    I2=double(hdf5read(image_2_filename,image_variable_name));
    
end;

% This determines the image size
image_size=zeros(1,3);
image_size(1)=size(I1,1);
image_size(2)=size(I1,2);
image_size(3)=size(I1,3);

% This iterates through the windows performing the cross-correlations







for window_index=1:prod(window_number);

    
    % This displays the calculation progress
    display_calculation_progress(window_index,1:size(window_min,1));
    
    % This extracts the correlation window from the first frame
    I1_Window=extract_correlation_window(I1,window_size,frame_1_window_min(window_index,:),frame_1_window_max(window_index,:),image_size);
    
    % This extracts the correlation window from the second frame
    I2_Window=extract_correlation_window(I2,window_size,frame_2_window_min(window_index,:),frame_2_window_max(window_index,:),image_size);
    
    % This subtracts the mean value from the windows if specified by the
    % user
    if zero_mean_windows;
        % This subtracts the mean value from the first image window
        I1_Window=I1_Window-mean(I1_Window(:));
        % This subtracts the mean value from the second image window
        I2_Window=I2_Window-mean(I2_Window(:));
    end;
    
   
    
    % This performs the cross-correlation measurement
    [Z,displacement_vector]=calculate_window_correlation(I1_Window,I2_Window,gaussian_filter,spectral_filter,window_size,correlation_method);
    
    % These are the velocities of the current correlation volume
    ii_velocity_current(window_index)=displacement_vector(1)+displacement_offset(window_index,1);
    jj_velocity_current(window_index)=displacement_vector(2)+displacement_offset(window_index,2);
    kk_velocity_current(window_index)=displacement_vector(3)+displacement_offset(window_index,3);

end;

% This overwrites the input velocity field coordinates with the current
% pass's coordinates
ii_position=ii_position_current;
jj_position=jj_position_current;
kk_position=kk_position_current;
% This overwrites the input velocity field vectors with the current
% pass's vectors
ii_velocity=ii_velocity_current/image_correlation_step;

size(ii_velocity_current)

jj_velocity=jj_velocity_current/image_correlation_step;
kk_velocity=kk_velocity_current/image_correlation_step;

%pyramid correlation location - originally

function [Z,displacement_vector]=calculate_window_correlation(I1,I2,gaussian_window_filter,spectral_filter,window_size,method);
% This performs the volumetric correlation between windows I1 and I2.  The
% gaussian filter is applied to the windows, the correlation is performed,
% and a sub-pixel fit is made to the correlation volume.

% This calculates the correlation volumes using either the SCC or RPC
% correlations
if strcmp(method,'SCC');

    
	% This performs the SCC correlation on the two windows
	[Z,C]=scc_correlation(I1,I2,gaussian_window_filter);
    
    
elseif strcmp(method,'RPC');

	% This performs the RPC correlation on the two windows
	C=rpc_correlation(I1,I2,gaussian_window_filter,spectral_filter);
    
end;



% This calculates the displacement vector from the correlation volumes
displacement_vector=subpixel_displacement_vector(C,window_size);




function [Z,C]=scc_correlation(I1,I2,gaussian_window_filter);
% This function performs the standard cross-correlation on the two windows 'I1' and 'I2'
% after masking the windows with a Gaussian function given by 'gaussian_window_filter'.

% This applies the Gaussian filter to the first window

I1=I1.*gaussian_window_filter;


% This applies the Gaussian filter to the second window
I2=I2.*gaussian_window_filter;

% This calculates the FFT of the first window
Y1=fftn(I1);

% This calculates the FFT of the second window
Y2=fftn(I2);

% This performs the cross-correlation in the spectral domain
Z= Y2.*conj(Y1);   % originally Y2.*conj(Y1)


% This performs the inverse FFT to return the correlation volume
C=ifftn(Z,'symmetric');




function C=rpc_correlation(I1,I2,gaussian_window_filter,spectral_filter);
% This function performs the robust phase correlation on the two windows 'I1' and 'I2'
% after masking the windows with a Gaussian function given by 'gaussian_window_filter' and
% using the RPC spectral filter defined by 'spectral_filter'.

% This applies the Gaussian filter to the first window
I1=I1.*gaussian_window_filter;

% This applies the Gaussian filter to the second window
I2=I2.*gaussian_window_filter;

% This calculates the FFT of the first window
Y1=fftn(I1);

% This calculates the FFT of the second window
Y2=fftn(I2);

% This performs the cross-correlation in the spectral domain
YC=Y2.*conj(Y1);

% This is the magnitude of the cross-correlation in the spectral domain
YC_Magnitude=sqrt(YC.*conj(YC));

% These are the indices of the non-zero magnitude components of the
% cross-correlation
non_zero_indices=(YC_Magnitude(:)~=0);

% This initializes the phase correlation array
R=YC;

% THis is the phase correlation in the spectral domain
R(non_zero_indices)=YC(non_zero_indices)./YC_Magnitude(non_zero_indices);

% This applies the spectral filter
R=R.*spectral_filter;

% This performs the inverse FFT to return the correlation volume
C=ifftn(R,'symmetric');

% This takes the absolute value of the correlation volume to ensure positive values
C=abs(C);



function displacement_vector=subpixel_displacement_vector(C,window_size);
% This function finds the subpixel fit to the correlation volume given by C.

% This is the linear index of the correlation maximum (this only takes the first maximum
% value if multiple are found - this is unlikely, but should be addressed)
max_linear_index=find(C==max(C(:)),1);

% This converts the linear index to a set of subscripted indices
[ii_max,jj_max,kk_max]=ind2sub(window_size,max_linear_index);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This performs the three-point Gaussian fitting in the ii-direction.                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% These are the three points to fit in the ii-direction
if (ii_max>1)&&(ii_max<window_size(1));
	% These are the three points to fit if the maximum is not on an edge
	C_ii_n1=C(ii_max-1,jj_max,kk_max);
	C_ii_00=C(ii_max,jj_max,kk_max);
	C_ii_p1=C(ii_max+1,jj_max,kk_max);
elseif (ii_max==1);
	% These are the three points to fit if the maximum is on the lower edge
	C_ii_n1=C(window_size(1),jj_max,kk_max);
	C_ii_00=C(ii_max,jj_max,kk_max);
	C_ii_p1=C(ii_max+1,jj_max,kk_max);
elseif (ii_max==window_size(1));
	% These are the three points to fit if the maximum is on the upper edge
	C_ii_n1=C(ii_max-1,jj_max,kk_max);
	C_ii_00=C(ii_max,jj_max,kk_max);
	C_ii_p1=C(1,jj_max,kk_max);
end;

% This is the numerator of the subpixel three-point Gaussian fit in the ii direction
ii_num_fit=log(C_ii_n1)-log(C_ii_p1);
% This is the denominator of the subpixel three-point Gaussian fit in the ii direction
ii_den_fit=2*log(C_ii_n1)-4*log(C_ii_00)+2*log(C_ii_p1);
% This is the subpixel three-point Gaussian fit in the ii direction
ii_fit=ii_max+ii_num_fit/ii_den_fit;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This performs the three-point Gaussian fitting in the jj-direction.                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% These are the three points to fit in the jj-direction
if (jj_max>1)&&(jj_max<window_size(2));
	% These are the three points to fit if the maximum is not on an edge
	C_jj_n1=C(ii_max,jj_max-1,kk_max);
	C_jj_00=C(ii_max,jj_max,kk_max);
	C_jj_p1=C(ii_max,jj_max+1,kk_max);
elseif (jj_max==1);
	% These are the three points to fit if the maximum is on the lower edge
	C_jj_n1=C(ii_max,window_size(2),kk_max);
	C_jj_00=C(ii_max,jj_max,kk_max);
	C_jj_p1=C(ii_max,jj_max+1,kk_max);
elseif (jj_max==window_size(2));
	% These are the three points to fit if the maximum is on the upper edge
	C_jj_n1=C(ii_max,jj_max-1,kk_max);
	C_jj_00=C(ii_max,jj_max,kk_max);
	C_jj_p1=C(ii_max,1,kk_max);
end;

% This is the numerator of the subpixel three-point Gaussian fit in the jj direction
jj_num_fit=log(C_jj_n1)-log(C_jj_p1);
% This is the denominator of the subpixel three-point Gaussian fit in the jj direction
jj_den_fit=2*log(C_jj_n1)-4*log(C_jj_00)+2*log(C_jj_p1);
% This is the subpixel three-point Gaussian fit in the jj direction
jj_fit=jj_max+jj_num_fit/jj_den_fit;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This performs the three-point Gaussian fitting in the kk-direction.                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This checks whether the array is 3D and if not, sets the sub-pixel fit
% value to 1 (which transforms to zero velocity)
if size(C,3)==1;
    
    % This sets the subpixel fit value in the 3rd dimension to one
    kk_fit=1;
    
else;

    % These are the three points to fit in the kk-direction
    if (kk_max>1)&&(kk_max<window_size(3));
        % These are the three points to fit if the maximum is not on an edge
        C_kk_n1=C(ii_max,jj_max,kk_max-1);
        C_kk_00=C(ii_max,jj_max,kk_max);
        C_kk_p1=C(ii_max,jj_max,kk_max+1);
    elseif (kk_max==1);
        % These are the three points to fit if the maximum is on the lower edge
        C_kk_n1=C(ii_max,jj_max,window_size(3));
        C_kk_00=C(ii_max,jj_max,kk_max);
        C_kk_p1=C(ii_max,jj_max,kk_max+1);
    elseif (kk_max==window_size(3));
        % These are the three points to fit if the maximum is on the upper edge
        C_kk_n1=C(ii_max,jj_max,kk_max-1);
        C_kk_00=C(ii_max,jj_max,kk_max);
        C_kk_p1=C(ii_max,jj_max,1);
    end;
    
    % This is the numerator of the subpixel three-point Gaussian fit in the kk direction
    kk_num_fit=log(C_kk_n1)-log(C_kk_p1);
    % This is the denominator of the subpixel three-point Gaussian fit in the kk direction
    kk_den_fit=2*log(C_kk_n1)-4*log(C_kk_00)+2*log(C_kk_p1);
    % This is the subpixel three-point Gaussian fit in the kk direction
    kk_fit=kk_max+kk_num_fit/kk_den_fit;
    
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This calculates the displacement vector from the sub-pixel fits.                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This accounts for the origin being located at the (1,1,1) index
ii_fit=ii_fit-1;
jj_fit=jj_fit-1;
kk_fit=kk_fit-1;

% This accounts for the periodicity of the correlation volume for the ii displacement
if ii_fit>(window_size(1)/2);
	ii_fit=ii_fit-window_size(1);
end;

% This accounts for the periodicity of the correlation volume for the jj displacement
if jj_fit>(window_size(2)/2);
	jj_fit=jj_fit-window_size(2);
end;

% This accounts for the periodicity of the correlation volume for the kk displacement
if kk_fit>(window_size(3)/2);
	kk_fit=kk_fit-window_size(3);
end;

% This is the displacement vector of the current window
displacement_vector=[ii_fit,jj_fit,kk_fit];



function I_Window=extract_correlation_window(I,window_size,window_min,window_max,image_size);
% This function extracts the current windows from the frame 'I' where 
% 'window_size' is the size of the correlation window, 'window_min' and 
% 'window_max' are the 1 x 3 vectors of the window minimum and maximum 
% indices, and 'image_size' is the the size of the frame 'I'.

% These are the coordinate domains to extract
ii_min=window_min(1);
ii_max=window_max(1);
jj_min=window_min(2);
jj_max=window_max(2);
kk_min=window_min(3);
kk_max=window_max(3);

% This checks if the ii indices are outside the domain of the image and sets them to the
% domain of the image if true
if ii_min<1;
	ii_min=1;
end;
if ii_max>image_size(1);
	ii_max=image_size(1);
end;
% This checks if the jj indices are outside the domain of the image and sets them to the
% domain of the image if true
if jj_min<1;
	jj_min=1;
end;
if jj_max>image_size(2);
	jj_max=image_size(2);
end;
% This checks if the kk indices are outside the domain of the image and sets them to the
% domain of the image if true
if kk_min<1;
	kk_min=1;
end;
if kk_max>image_size(3);
	kk_max=image_size(3);
end;

% This extracts the current window from the frame
I_Window=I(ii_min:ii_max,jj_min:jj_max,kk_min:kk_max);

% This initilizes the image size vector
I_Window_Size=zeros(1,3);
% This is the size of the extracted image
I_Window_Size(1)=size(I_Window,1);
I_Window_Size(2)=size(I_Window,2);
I_Window_Size(3)=size(I_Window,3);

% If the extracted image is smaller than the full window size (due to the window being
% on the edge of the full image) this pads the image with zeros
if any(I_Window_Size~=window_size);
	% This initializes an array of zeros the size of the full window
	I_Window_Temp=zeros(window_size);
	% These are the indices to insert the image into the zero padded image (the
	% indices are centered in the window so that information is not lost during the
	% masking process)
	window_insert_min=floor((window_size-I_Window_Size)/2)+1;   	
	window_insert_max=window_insert_min+I_Window_Size-1;
	% This copies the current image from the frame into the window
	I_Window_Temp(window_insert_min(1):window_insert_max(1),window_insert_min(2):window_insert_max(2),window_insert_min(3):window_insert_max(3))=I_Window;
	% This renames the variable so that the zero-padded image is used
	I_Window=I_Window_Temp;
end;




function p=determine_gaussian_size_2(window_size,resolution_size);
% This function is used to determine the diameter of a Gaussian function
% used to mask PIV windows to prevent aliasing during the Fourier
% transform.
%
% The input argument 'window_size' refers to the full size of
% the PIV window prior to performing the masking operation and may be a
% scaler or a vector.  The input argument 'resolution_size' refers to the
% effective resolution that the user wants to attain from the masked PIV
% window, ie the window may be 64 x 64 in size (window_size=[64,64]) but
% the user wants the effective resolution after masking to be equal to 32 x
% 32 (resolution_size=[32,32]).
%
% The output argument 'p' is used to contruct the Gaussian masking function
% of the following form
%
%  G(x) = exp(-((p*x)^2)/2)                                             (1)
%
% where the independent variable lies in the domain
%
%  -(window_size-1)/2 <= x <= (window_size-1)/2                         (2)
%
% The variable 'p' is calculated by solving the integral equation
%
%  resolution_size = int(G(x),-(window_size-1)/2,(window_size-1)/2)     (3)
%
% or the equivalent expression
%
%  resolution_size = sqrt(2*pi)*erf(window_size*p/(2*sqrt(2))/p         (4)
%
% given by simplifying the integral.  These equations are taken from
% "Assessment of advanced windowing techniques for digital particle image
% velocimetry (DPIV)" Eckstein 2009.
%
% This function is written to replace the 'findwidth' function in prana.
% The equation is solved by re-writing equation (4) in the form
% 
%  0 = f(x) = erf(b*x)/x-a                                              (5)
% 
% and solving for the independent variable using Halley's method.
%
% Written By: Rod La Foy
% Written On: 6 June 2013

% These are the ratios of the resolution size to window size
size_ratio=resolution_size./window_size;

% This initializes the output vector p
p=zeros(size(window_size));

% This pre-computes the square root of pi for later use
pi_sqrt=sqrt(pi);

% This is the tolerance with which to find the root
tolerance=1e-4;

% This iterates through the number of dimensions to solve for (probably
% either 2 or 3)
for dimension_number=1:length(window_size);

	% This checks if the size ratio is equal to one or greater
	if size_ratio(dimension_number)>=1;
		
		% This sets the Gaussian width variable to 0
		p(dimension_number)=0;
		
	else;
    
		% This is the first constant term in the equation
		a=resolution_size(dimension_number)/sqrt(2*pi);
		% This is the second constant term in the equation
		b=window_size(dimension_number)/(2*sqrt(2));

		% This is the initial guess of value of the Gaussian width term
		p_temp=1;
	
		% This iterates to the root
		while true;
		
			% This is the error function of the b*p term
			bp_erf=erf(b*p_temp);
			% This is the exponential term of the (b*p)^2 term
			bp_exp=exp((b*p_temp)^2);
		
			% This is the numerator of the difference term
			dp_numerator=a*p_temp-bp_erf;
		
			% This initializes the denominator term
			dp_denominator=0;
			% This adds the first term to the denominator
			dp_denominator=dp_denominator+a;
			% This adds the second term to the denominator
			dp_denominator=dp_denominator+2*b/(pi_sqrt*bp_exp);
			% This adds the third term to the denominator
			dp_denominator=dp_denominator+(2*b^3*p_temp^2*dp_numerator)/(2*b*p_temp-bp_exp*pi_sqrt*bp_erf);
		
			% This is the expected difference to the root
			dp=dp_numerator/dp_denominator;
		
			% This is the new root approximation
			p_temp=p_temp-dp;
		
			% This breaks if the error is less then the tolerance
			if abs(erf(b*p_temp)/p_temp-a)<tolerance;
				% This breaks the loop
				break;
			end;
		
		end;
	
		% This is the Gaussian function diameter term
		p(dimension_number)=p_temp;
		
	end;
    
end;



function gaussian_window_filter=create_gaussian_filter(window_size,gaussian_width);
% This function creates a Gaussian filter to apply in the spatial domain to
% the correlation windows to remove the effects of aliasing

% This is the standard deviation of the filter
sigma=1./gaussian_width;

% These are the vectors of coordinates used to create the Gaussian function
ii_index_vector=linspace(-(window_size(1)-1)/2,(window_size(1)-1)/2,window_size(1));
jj_index_vector=linspace(-(window_size(2)-1)/2,(window_size(2)-1)/2,window_size(2));
kk_index_vector=linspace(-(window_size(3)-1)/2,(window_size(3)-1)/2,window_size(3));

% These are the full array of coordinates
[ii_index_array,jj_index_array,kk_index_array]=ndgrid(ii_index_vector,jj_index_vector,kk_index_vector);

% This creates the Gaussian function
gaussian_window_filter=exp(-(ii_index_array.^2)/(2*sigma(1)^2)-(jj_index_array.^2)/(2*sigma(2)^2)-(kk_index_array.^2)/(2*sigma(3)^2));



function [frame_1_window_min,frame_1_window_max,frame_2_window_min,frame_2_window_max,displacement_offset]=dwo_window_domains(window_min,window_max,window_center,image_correlation_step,ii_coordinates,jj_coordinates,kk_coordinates,ii_displacement,jj_displacement,kk_displacement);
% This function adds the discrete window offset to the window positions
% based upon the known velocity field.  The non-offset window positions are
% given by the variables
%
%  'window_min'
%  'window_max'
%  'window_center'
%
% where these variables are N x 3 arrays where N is the total number of
% windows and column vectors correspond to the 1st, 2nd, and 3rd dimension
% coordinates.  The known velocity field is given at the coordinates
%
%  'ii_coordinates'
%  'jj_coordinates'
%  'kk_coordinates'
%
% which are I x J x K arrays where I is the number of known vector
% coordinates in the 1st dimension, J is the number of known vector
% coordinates in the 2nd dimension, and K is the number of known vector
% coordinates in teh 3rd dimension.  The known velocity field values are
% given by the arrays
%
%  'ii_displacement'
%  'jj_displacement'
%  'kk_displacement'
%
% which are also I x J x K in size.  The output arrays
%
%  'frame_1_window_min'
%  'frame_1_window_max'
%  'frame_2_window_min'
%  'frame_2_window_max'
%
% are N x 3 vectors in size and given the minimum and maximum range of the
% first and second frames respectively.  The output array
%
%  'displacement_offset'
%
% is also N x 3 in size and corresponds to the offset that must be added
% back to the measured velocity field due to the discrete window offset.

% This scales the velocity estimates by the number of number of frames
% between the correlations
ii_displacement=image_correlation_step*ii_displacement;
jj_displacement=image_correlation_step*jj_displacement;
kk_displacement=image_correlation_step*kk_displacement;

% This checks whether the array of velocity values is 2D or 3D and chooses
% the appropriate interpolotion function
if size(ii_coordinates,3)==1;
    
    % If the vector field is greater then 3 x 3, then a cubic interpolation
    % is performed, otherwise a linear interpolation is performed
    if (size(ii_coordinates,1)>2)&&(size(ii_coordinates,2)>2);
        
        % This is the estimated ii-displacement at the center of the windows
        estimated_window_dii=interp2(ii_coordinates,jj_coordinates,ii_displacement,window_center(:,1),window_center(:,2),'cubic',0);
        % This is the estimated jj-displacement at the center of the windows
        estimated_window_djj=interp2(ii_coordinates,jj_coordinates,jj_displacement,window_center(:,1),window_center(:,2),'cubic',0);
        % This is the estimated kk-displacement at the center of the windows
        estimated_window_dkk=interp2(ii_coordinates,jj_coordinates,kk_displacement,window_center(:,1),window_center(:,2),'cubic',0);
        
    else;
        
        % This is the estimated ii-displacement at the center of the windows
        estimated_window_dii=interp2(ii_coordinates,jj_coordinates,ii_displacement,window_center(:,1),window_center(:,2),'linear',0);
        % This is the estimated jj-displacement at the center of the windows
        estimated_window_djj=interp2(ii_coordinates,jj_coordinates,jj_displacement,window_center(:,1),window_center(:,2),'linear',0);
        % This is the estimated kk-displacement at the center of the windows
        estimated_window_dkk=interp2(ii_coordinates,jj_coordinates,kk_displacement,window_center(:,1),window_center(:,2),'linear',0);
        
    end;
    
else;
    
    % If the vector field is greater then 3 x 3, then a cubic interpolation
    % is performed, otherwise a linear interpolation is performed
    if (size(ii_coordinates,1)>2)&&(size(ii_coordinates,2)>2)&&(size(ii_coordinates,3)>2);
        
        % This is the estimated ii-displacement at the center of the windows
        estimated_window_dii=interp3(ii_coordinates,jj_coordinates,kk_coordinates,ii_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0);
        % This is the estimated jj-displacement at the center of the windows
        estimated_window_djj=interp3(ii_coordinates,jj_coordinates,kk_coordinates,jj_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0);
        % This is the estimated kk-displacement at the center of the windows
        estimated_window_dkk=interp3(ii_coordinates,jj_coordinates,kk_coordinates,kk_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0);
        
    else;
        
        % This is the estimated ii-displacement at the center of the windows
        estimated_window_dii=interp3(ii_coordinates,jj_coordinates,kk_coordinates,ii_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0);
        % This is the estimated jj-displacement at the center of the windows
        estimated_window_djj=interp3(ii_coordinates,jj_coordinates,kk_coordinates,jj_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0);
        % This is the estimated kk-displacement at the center of the windows
        estimated_window_dkk=interp3(ii_coordinates,jj_coordinates,kk_coordinates,kk_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0);
        
    end;

end;
    
% This is half of the measured velocity at the window coordinates
rounded_half_estimated_displacement=round([estimated_window_dii,estimated_window_djj,estimated_window_dkk]/2);

% This calculates the first frame minimum domain values
frame_1_window_min=window_min-floor(rounded_half_estimated_displacement);
% This calculates the first frame maximum domain values
frame_1_window_max=window_max-floor(rounded_half_estimated_displacement);

% This calculates the second frame minimum domain values
frame_2_window_min=window_min+ceil(rounded_half_estimated_displacement);
% This calculates the second frame maximum domain values
frame_2_window_max=window_max+ceil(rounded_half_estimated_displacement);

% This is the offset that must be added back onto the measured velocity
% field of the current pass to account for the dicrete window offset
displacement_offset=ceil(rounded_half_estimated_displacement)+floor(rounded_half_estimated_displacement);



function [window_min,window_max,window_center,window_number]=calculate_window_domains(image_size,window_resolution,window_size,grid_spacing);
% This function calculates the minimum and maximum values of the window
% indices to extract from the image.

% These are the number of windows in each dimension
window_number=ceil((image_size-window_resolution)./grid_spacing); %ORIGINALLY window_number=ceil((image_size-window_resolution)./grid_spacing)+1

% This initializes the window indices
window_min=zeros(prod(window_number),3);
window_max=zeros(prod(window_number),3);

% This initializes the window centers
window_center=zeros(prod(window_number),3);

% This is the full resolution covered by the PIV windows (which will be equal to or larger
% than the resolution of the full image)
full_resolution=grid_spacing.*(window_number-1)+window_resolution;





% This is the distance in indices below 1 to start the windows
negative_index_distance=floor((full_resolution-image_size)/2);

% This is a counting variable for indexing the window domains
count=0;
	
% This calculates the window domains
for kk=1:window_number(3);
	for jj=1:window_number(2);
		for ii=1:window_number(1);
		
			% This increments the counting variable
			count=count+1;
			
			% These are the minimum window indices
			window_min(count,:)=([ii,jj,kk]-1).*grid_spacing-ceil((window_size-window_resolution)/2)-negative_index_distance+1;
			% These are the maximum window indices
			window_max(count,:)=window_min(count,:)+window_size-1;

			% These are the window center indices
            window_center(count,:)=(window_min(count,:)+window_max(count,:))/2 - 0.5;

            window_min(count,:)
            window_max(count,:)

		end;
	end;
end;



function spectral_filter=create_spectral_filter(window_size,particle_diameter);
% This function creates the RPC spectral energy filter with a size of 'window_size' using
% a particle diameter given by 'particle_diameter' where 'bit_depth' is the bit depth of
% the three-dimensional image to be processed.

% These are the wavenumber vectors
k_ii_vector=-pi:2*pi/window_size(1):pi-2*pi/window_size(1);
k_jj_vector=-pi:2*pi/window_size(2):pi-2*pi/window_size(2);
k_kk_vector=-pi:2*pi/window_size(3):pi-2*pi/window_size(3);

% These are the wavenumber coordinates
[k_ii,k_jj,k_kk]=ndgrid(k_ii_vector,k_jj_vector,k_kk_vector);

% This is the squared radius of the wavenumber
rho=k_ii.^2+k_jj.^2+k_kk.^2;

% None of the coefficients actually matter . . . so RPC filter is just a low pass filter . . .
spectral_filter=exp(-(rho*particle_diameter^2)/16);



function display_calculation_progress(current_value,value_vector);
% This function displays a progress bar showing the percent complete of the
% currently running calculation where the calculation is being iterated
% over the vector 'value_vector' and the current iteration's value is equal
% to 'current_value'.  For example in the for loop
%
%  for ii=1:10;
%       commands . . .
%  end;
%
% the arguments would be equal to
%
%  current_value=ii;
%  value_vector=1:10;
%
% This convention holds for non-integer or non-monotonic vectors, although
% to work correctly all values of value_vector must be unique.

% This is the number of characters to display in the pogress bar (denoted
% by either '#' or '-' characters)
progress_bar_character_number=50;

% This is the index of the value_vector that corresponds to the
% current_value
[null_variable,value_index]=min(abs(current_value-value_vector));

% This is the percentage of the calculation that is completed
current_progress_decimal=value_index/length(value_vector);
% This creates the character string showing the numerical and graphical
% progress for the current iteration
current_text_string=generate_progress_string(current_progress_decimal,progress_bar_character_number);

% If this is the first iteration, then a new line is added, the progress
% bar is displayed, and the function is exited
if value_index==1;
    % This displays the portion of the string showing the percentage of 
    % the calculation that is complete that is new to this iteration
    fprintf(current_text_string);
    % This ends the function
    return;
end;

% This is the percentage of the calculation that was completed during the
% last iteration
previous_progress_decimal=(value_index-1)/length(value_vector);
% This creates the character string showing the numerical and graphical
% progress for the previous iteration
previous_text_string=generate_progress_string(previous_progress_decimal,progress_bar_character_number);

% This compares the current progress string with the previous string and if
% they are the same, then the function exits.  If they are different, then
% only the text after the difference is displayed.
if strcmp(current_text_string,previous_text_string);
    
    % If this is the last time that the progress bar will be displayed, this
    % prints a new line
    if value_index==length(value_vector);
        % This prints a new line after the progress bar
        fprintf('\n');
    end;
    
    % This exits the function without changing the progress bar
    return;
    
else;
    % This is the total number of charcters to be displayed
    string_character_number=length(current_text_string);
    
    % This is the index into the string where the strings first differ
    first_difference_index=find(not(current_text_string==previous_text_string),1,'first');
    
    % These are the locations of the double percent signs '%%'
    double_percent_indices=strfind(current_text_string,'%%');
    % This removes the double percent indices that are before the first
    % difference index
    double_percent_indices(double_percent_indices<first_difference_index)=[];

    % This is the number of characters of the previous line to delete
    delete_character_number=string_character_number-first_difference_index+1-length(double_percent_indices);
    % If this is the first iteration, then no characters are deleted
    if value_index==1;
        % This sets the number of characters to be deleted to zero
        delete_character_number=0;
        % This sets the first difference character to one
        first_difference_index=1;
    end;
    
    % This deletes the previously displayed characters back to the first
    % differing character (by displaying the 'backspace' character)
    fprintf(1,repmat('\b',1,delete_character_number));
    
    % This displays the portion of the string showing the percentage of 
    % the calculation that is complete that is new to this iteration
    fprintf(current_text_string(first_difference_index:end));
    
    % If this is the last time that the progress bar will be displayed, this
    % prints a new line
    if value_index==length(value_vector);
        % This prints a new line after the progress bar
        fprintf('\n');
    end;
    
end;



function text_string=generate_progress_string(progress_decimal,progress_bar_character_number);
% This function generates the progress bar text that will be displayed for
% the current percentage of the calculation that is completed.

% This is a string giving a numerical value of the percentage of the
% calculation completed
numerical_progress_string=sprintf('% 4.0f',round(100*progress_decimal));

% This is the prefix to the progress bar showing the numerical value of the
% percent complete
string_prefix=['Calculation',numerical_progress_string,'%% Complete:     0%% ['];

% This is the suffix to the progress bar
string_suffix='] 100%%';

% This is the number of '#' signs to display corresponding to the graphical
% representation of the percent of the calculation that is complete
completed_character_number=round(progress_decimal*progress_bar_character_number);
% This is the number of '-' signs to display corresponding to the graphical
% representation of the percent of the calculation that remains
remaining_character_number=progress_bar_character_number-completed_character_number;

% This creates the string of characters representing the graphical
% percentage of the calculation that is complete
progress_bar_string=[repmat('#',1,completed_character_number),repmat('-',1,remaining_character_number)];

% This is the total text string to be displayed
text_string=[string_prefix,progress_bar_string,string_suffix];