diff --git a/PIV_2D_wdef.m b/PIV_2D_wdef.m
index fedce4c24696491ee07846b641514c316ac12d6a..1dfd7cffa70f17b131cfd5f82c571f2886989766 100644
--- a/PIV_2D_wdef.m
+++ b/PIV_2D_wdef.m
@@ -152,6 +152,11 @@ function piv_result = PIV_2D_wdef(A, B, mask, piv_setup)
% window size: we need to size this filter relative to the cutoff frequency
ksize_filt = round(wsize_old ./ win_spacing_old) + 1;
sd = sqrt(prod(ksize_filt))/3*0.65;
+ % median filter to avoid errors due to spurious vectors in
+ % result
+ %delta_ab_old(:,:,1) = medfilt2(delta_ab_old(:,:,1), ksize_filt);
+ %delta_ab_old(:,:,2) = medfilt2(delta_ab_old(:,:,2), ksize_filt);
+ % smoothing filter
G_smooth = fspecial('gaussian', ksize_filt, sd);
delta_ab_old = convn(delta_ab_old, G_smooth, 'same');
@@ -284,9 +289,19 @@ function piv_result = PIV_2D_wdef(A, B, mask, piv_setup)
% mark bad vectors as NaN
ux_mat(nan_mask) = nan;
uy_mat(nan_mask) = nan;
- % inpaint
- ux_mat = inpaint_nans(ux_mat,1);
- uy_mat = inpaint_nans(uy_mat,1);
+
+ if true % pass == 1
+ % inpaint NaN values on first pass
+ ux_mat = inpaint_nans(ux_mat,1);
+ uy_mat = inpaint_nans(uy_mat,1);
+ else
+ % subsequent passes just us previous predictor field
+ ux_pred = delta_ab_pred(:,:,1);
+ uy_pred = delta_ab_pred(:,:,2);
+ ux_mat(nan_mask)= ux_pred(nan_mask);
+ uy_mat(nan_mask)= uy_pred(nan_mask);
+ end
+
peak_choice(nan_mask)= 0;
Q(nan_mask) = 0;
peak_mag(nan_mask) = 0;
diff --git a/calibration/chessboard_calibration.m b/calibration/chessboard_calibration.m
new file mode 100644
index 0000000000000000000000000000000000000000..61e3f5739ffd8a270633b1d31df66cb706aaccf1
--- /dev/null
+++ b/calibration/chessboard_calibration.m
@@ -0,0 +1,250 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% dotboard_calibration.m
+% find markers in images of calibration plates
+% then perform a 2D calibration to dewarp the calibration plate image
+
+addpath([pwd '/../minmax/']);
+addpath([pwd '/../pinhole/']);
+addpath([pwd '/../misc/']);
+
+%% input parameters
+
+grid_dx = -1; % real units, mm
+grid_dy = +1;
+invert_image = 0; % invert image for black dots on white
+blur_image = 0; % gaussian blur image
+LoG_image = 0; % apply laplacian of Gaussian filter
+checkerboard_image = 0; % apply checkerboard image filter
+median_filter = 1;
+median_filter_size = [5 5];
+normalise_image = 1;
+corner_metric = 0.1; % threshold for identifying corners
+dot_half_window = [64 64];
+
+
+output_dir = 'E:\mixer\lpt\calibration\21-06-14\';
+input_file_list = {'E:\mixer\lpt\calibration\21-06-14\C15794_im02.tif'};
+plate_origin = [0 0 0]';
+camera_index = 1;
+camera_id = 1;
+calib_idx_offset = 0;
+origin_pos = [0 0];
+
+image_clim = [0 255];
+
+%% generate calibration for each image
+
+nCalibrations = length(input_file_list);
+
+calib_idx = 1;
+
+while calib_idx <= nCalibrations
+ %% update user of status
+ input_file = input_file_list{calib_idx};
+ output_file = sprintf('%s/CAM%d-%02d.mat', output_dir, camera_index, calib_idx+calib_idx_offset);
+ fprintf('Image %02d of %02d:\n', calib_idx, nCalibrations);
+
+ % check if calibration already exists
+ if exist(output_file, 'file')
+ sInput = input('Calibration already exists - overwrite? y to overwrite: ', 's');
+ if ~strcmpi(sInput, 'y')
+ calib_idx = calib_idx + 1;
+ continue;
+ end
+ end
+ % check if image input exists
+ if ~exist(input_file, 'file')
+ fprintf(' - could not find image file\n');
+ end
+
+ %% load image
+ im = readimg(input_file);
+ im = double((im{1}'));
+ Nx = size(im, 1);
+ Ny = size(im, 2);
+ i_axis = 0 : Nx-1;
+ j_axis = 0 : Ny-1;
+
+ %% create filtered image for processing
+ im_filt = im;
+ if invert_image
+ im_filt = max(im(:))-im;
+ end
+ if blur_image
+ G = fspecial('gaussian', [5 5]);
+ im_filt = filter2(G, im_filt);
+ end
+ if LoG_image
+ scale = norm(dot_half_window)*0.2;
+ G = fspecial('log', 2*dot_half_window, scale);
+ im_filt = scale * conv2(G, im_filt);
+ end
+ if median_filter
+ im_filt = medfilt2(im_filt, median_filter_size);
+ end
+ if checkerboard_image
+ G = repelem(eye(2), dot_half_window(1), dot_half_window(2));
+ im_filt = conv2(G, im_filt);
+ end
+ if normalise_image
+ [minv, maxv] = minmaxfiltnd(single(im_filt), dot_half_window*2+1);
+ im_filt = (im_filt - minv) ./ max(16, maxv-minv);
+ end
+
+ %% display image
+ figure(1);
+ clf;
+ hold off;
+ ax=gca;
+ h_filt = imagesc(im_filt');
+ axis equal tight xy;
+ colormap gray;
+ hold on;
+ colorbar;
+
+ %% crop image
+ fprintf('Select the bottom left square of the board:\n');
+ im_bl = ginput(1);
+ plot(im_bl(1), im_bl(2), 'ro');
+ fprintf('Select the bottom right square of the board:\n');
+ im_br = ginput(1);
+ plot(im_br(1), im_br(2), 'ro');
+ fprintf('Select the top right square of the board:\n');
+ im_tr = ginput(1);
+ plot(im_tr(1), im_tr(2), 'ro');
+ fprintf('Select the top left square of the board:\n');
+ im_tl = ginput(1);
+ plot(im_tl(1), im_tl(2), 'ro');
+ % identify corners
+ im_corners = [im_bl; im_tr];
+ im_corners = sort(round(im_corners), 1);
+
+ %% coordinate directions
+ im_ex = im_br + im_tr - im_bl - im_tl;
+ im_ey = im_tr + im_tl - im_br - im_bl;
+ im_ex = im_ex' / norm(im_ex);
+ im_ey = im_ey' / norm(im_ey);
+
+ %% display cropped and filtered image
+ plot(ax(1), im_corners([1 2 2 1 1], 1), im_corners([1 1 2 2 1], 2), 'r-');
+ rx = i_axis > im_corners(1,1) & i_axis < im_corners(2,1);
+ ry = j_axis > im_corners(1,2) & j_axis < im_corners(2,2);
+ % crop filtered image
+ im_crop = zeros(Nx,Ny);
+ im_crop(rx,ry) = im_filt(rx,ry);
+ im_filt = im_crop;
+ set(h_filt, 'cdata', im_filt');
+ drawnow;
+
+ %% extract markers using chessboard algorithm
+ h_markers = [];
+ h_text = [];
+ while ~isnan(corner_metric)
+ %% search for chessboard corners
+ fprintf('searching for checkerboard corners with m=%0.3f ...\n', corner_metric);
+ [ij_markers, board_size] = detectCheckerboardPoints(im_filt', 'MinCornerMetric', corner_metric);
+ ij_markers = ij_markers'; % 2 x n_markers
+
+ %% update figure or try again
+ if isempty(ij_markers)
+ fprintf('no valid markers found.\n');
+ fprintf('try again with a different threshold\n\n');
+ else
+ %% display
+ figure(1);
+ hold on;
+ delete(h_markers);
+ delete(h_text);
+ h_markers = plot(ax(1), ij_markers(1,:), ij_markers(2,:), 'g+');
+ end
+
+ %% prompt for new threshold
+ ret = input('enter a new corner metric or nan to accept: ', 's');
+ corner_metric = str2double(ret);
+ end
+
+ %% get coordinates of markers
+ % create grid
+ x_grid = (1 : board_size(2)-1) * grid_dx;
+ y_grid = (1 : board_size(1)-1) * grid_dy;
+ x_grid = x_grid - round(mean(x_grid));
+ y_grid = y_grid - round(mean(y_grid));
+ [xm, ym] = meshgrid(x_grid, y_grid);
+ xy_markers = [xm(:) ym(:)]';
+ n_markers = prod(board_size-1);
+
+ % annotate
+ h_text = [];
+ delete(findobj(gcf,'type','text'));
+ for k = [1 board_size(1)-1 n_markers n_markers-board_size(1)+2]
+ str = sprintf('(%0.1f,%0.1f)', xy_markers(:,k));
+ h_text(k) = text(ij_markers(1,k), ij_markers(2,k), str, 'color', 'b', ...
+ 'fontsize', 18, ...
+ 'horizontalalignment', 'center', 'verticalalignment', 'middle');
+ end
+
+ %% create calibration and dewarp
+ xy_tilde = [xy_markers; ones(1,n_markers)];
+ P2D = dlt33([ij_markers;ones(1,n_markers)], xy_tilde);
+ ij_tilde = P2D * xy_tilde;
+ lambda = [1; 1] * ij_tilde(3,:);
+ ij_fit = ij_tilde(1:2,:) ./ lambda;
+
+ ij_error = ij_markers - ij_fit;
+ px_error = sqrt(ij_error(1,:).^2 + ij_error(2,:).^2);
+ rms_error = sqrt(mean(px_error.^2));
+
+ fprintf('RMS fit error: %0.2f px\n', rms_error);
+
+ %% work out x and y axes
+ x_lim = sort(x_grid([1 end]));
+ y_lim = sort(y_grid([1 end]));
+ x_axis = linspace(x_lim(1), x_lim(2), Nx);
+ y_axis = linspace(y_lim(1), y_lim(2), Ny);
+ % calculate dewarping
+ [xmat, ymat] = ndgrid(x_axis, y_axis);
+ ij_tilde_dewarped = P2D * [xmat(:) ymat(:) ones(Nx*Ny,1)]';
+ lambda = ij_tilde_dewarped(3,:);
+ imat = reshape(ij_tilde_dewarped(1,:) ./ lambda, [Nx Ny]);
+ jmat = reshape(ij_tilde_dewarped(2,:) ./ lambda, [Nx Ny]);
+ im_xy = interp2(i_axis, j_axis, im', imat, jmat, 'linear');
+
+ %% show dewarped image
+ figure(2);
+ hold off;
+ imagesc(x_axis, y_axis, im_xy');
+ xlim(x_lim);
+ ylim(y_lim);
+ colormap gray;
+ axis xy equal tight;
+ hold on;
+ %set(gca,'Clim', image_clim);
+ for x = (round(x_lim(1)/abs(grid_dx)) : round(x_lim(end)/abs(grid_dx))) * abs(grid_dx)
+ plot([x x], [y_axis(1) y_axis(end)], 'r-');
+ end
+ for y = (round(y_lim(1)/abs(grid_dy)) : round(y_lim(end)/abs(grid_dy))) * abs(grid_dy)
+ plot([x_axis(1) x_axis(end)], [y y], 'r-');
+ end
+
+ %% query user to save calibration or repeat
+ sInput = input('Would you like to repeat the calibration? y to repeat: ', 's');
+ if strcmpi(sInput, 'y')
+ continue; % continue without incrementing the counter
+ end
+
+ %% save output
+ fprintf('Saving ... ');
+ fstruct = struct('source_image', im, ...
+ 'dewarped_image', im_xy, ...
+ 'plate_pos', plate_origin(:, calib_idx), ...
+ 'P', P2D, ...
+ 'i_axis', i_axis, 'j_axis', j_axis, ...
+ 'x_axis', x_axis, 'y_axis', y_axis, 'x_lim', x_lim, 'y_lim', y_lim, ...
+ 'grid_dx', grid_dx, 'grid_dy', grid_dy, ...
+ 'ij_markers', ij_markers, 'xy_markers', xy_markers, ...
+ 'camera_index', camera_index, ...
+ 'id', camera_id);
+ save(output_file, '-struct', 'fstruct');
+ fprintf('done\n');
+ calib_idx = calib_idx + 1;
+end
\ No newline at end of file
diff --git a/calibration/dotboard_calibration.m b/calibration/dotboard_calibration.m
new file mode 100644
index 0000000000000000000000000000000000000000..b800705b8c138f950a4cbef668ff3a2c9a6d0d00
--- /dev/null
+++ b/calibration/dotboard_calibration.m
@@ -0,0 +1,222 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% dotboard_calibration.m
+% find markers in images of calibration plates
+% then perform a 2D calibration to dewarp the calibration plate image
+
+addpath([pwd '/../../mex/minmax/']);
+addpath([pwd '/pinhole/']);
+
+%% input parameters
+
+grid_dx = 20; % real units, mm
+grid_dy = 20;
+invert_image = 1; % invert image for black dots on white
+blur_image = 0; % gaussian blur image
+LoG_image = 0; % apply laplacian of Gaussian filter
+normalise_image = 0;
+
+
+output_dir = 'M:\mixer\calibration\20-11-17';
+input_file_list = { 'M:\mixer\calibration\20-11-17\B00001.tif'};
+plate_origin = [0 0 0]';
+camera_index = 1;
+camera_id = 1;
+dot_half_window = [32 32];
+calib_idx_offset = 0;
+origin_pos = [0 0];
+
+image_clim = [0 255];
+
+%% generate calibration for each image
+
+nCalibrations = length(input_file_list);
+
+calib_idx = 1;
+
+while calib_idx <= nCalibrations
+ %% update user of status
+ input_file = input_file_list{calib_idx};
+ output_file = sprintf('%s/CAM%d-%02d.mat', output_dir, camera_index, calib_idx+calib_idx_offset);
+ fprintf('Image %02d of %02d:\n', calib_idx, nCalibrations);
+
+ % check if calibration already exists
+ if exist(output_file, 'file')
+ sInput = input('Calibration already exists - overwrite? y to overwrite: ', 's');
+ if ~strcmpi(sInput, 'y')
+ calib_idx = calib_idx + 1;
+ continue;
+ end
+ end
+ % check if image input exists
+ if ~exist(input_file, 'file')
+ fprintf(' - could not find image file\n');
+ end
+
+ %% load image, request position
+ im = readimg(input_file);
+ %im = double(im{2}');
+ im = rot90(im{2},-1);
+ Nx = size(im, 1);
+ Ny = size(im, 2);
+
+ %% create filtered image for processing
+ im_filt = im;
+ if invert_image
+ im_filt = max(im(:))-im;
+ end
+ if blur_image
+ G = fspecial('gaussian', [5 5]);
+ im_filt = filter2(G, im_filt);
+ end
+ if LoG_image
+ scale = norm(dot_half_window)*0.2;
+ G = fspecial('log', 2*dot_half_window, scale);
+ im_filt = scale * conv2(G, im_filt);
+ end
+ if normalise_image
+ [minv, maxv] = minmaxfiltnd(single(im_filt), dot_half_window*2+1);
+ im_filt = (im_filt - minv) ./ max(16, maxv-minv);
+ end
+
+ %% display image
+ figure(1);
+ hold off;
+ imagesc(im_filt');
+ axis equal tight xy;
+ colormap gray;
+ %set(gca,'Clim', image_clim);
+ hold on;
+ colorbar;
+
+ %% select region of interest
+ fprintf('Please select region of interest');
+ [i_coord, j_coord] = ginput(2);
+ xlim(i_coord);
+ ylim(j_coord);
+
+ %% get six marker positions
+ % center, first right, far right, far corner, far above, above
+ bDone = false;
+ while ~bDone
+ %% get marker position from user
+ fprintf('Please identify the six markers\n');
+ fprintf('Press any key to abort and try again\n');
+
+ i_marker = nan(1,6);
+ j_marker = nan(1,6);
+ b_aborted = false;
+ h_plot = [];
+ for mark = 1 : 6
+ %% get marker position from user input
+ [i_marker(mark) j_marker(mark), btn] = ginput(1);
+ if btn > 3
+ b_aborted= true;
+ break;
+ end
+ %% get barycentre of marker, to be more accurate
+ dot_img_i = round(i_marker(mark)) + (-dot_half_window(1) : dot_half_window(1));
+ dot_img_j = round(j_marker(mark)) + (-dot_half_window(2) : dot_half_window(2));
+ [imat,jmat] = ndgrid(dot_img_i, dot_img_j);
+ dot_img = im_filt(dot_img_i, dot_img_j);
+ i_marker(mark) = sum(dot_img(:) .* imat(:)) / sum(dot_img(:));
+ j_marker(mark) = sum(dot_img(:) .* jmat(:)) / sum(dot_img(:));
+
+ %% show user what he's plotting
+ delete(h_plot);
+ h_plot = plot(i_marker, j_marker, 'ro-');
+ for it = 1 : mark
+ plot( i_marker(it)+[-1 -1 +1 +1 -1]*dot_half_window(1), ...
+ j_marker(it)+[-1 +1 +1 -1 -1]*dot_half_window(2), 'g-');
+ end
+ end
+ if b_aborted
+ delete(h_plot);
+ continue;
+ end
+ bDone = true;
+ end
+
+ %% find marker positions
+ [ij_markers, xy_markers]= find_markers(im_filt, [i_marker;j_marker], dot_half_window);
+ xy_markers(1,:) = xy_markers(1,:) * grid_dx + origin_pos(1);
+ xy_markers(2,:) = xy_markers(2,:) * grid_dy + origin_pos(2);
+ n_markers = size(xy_markers,2);
+ xy_tilde = [xy_markers; ones(1,n_markers)];
+
+ %% create calibration and dewarp
+ P2D = dlt33([ij_markers;ones(1,n_markers)], xy_tilde);
+ ij_tilde = P2D * xy_tilde;
+ lambda = [1; 1] * ij_tilde(3,:);
+ ij_fit = ij_tilde(1:2,:) ./ lambda;
+
+ ij_error = ij_markers - ij_fit;
+ px_error = sqrt(ij_error(1,:).^2 + ij_error(2,:).^2);
+ rms_error = sqrt(mean(px_error.^2));
+
+ fprintf('RMS fit error: %0.2f px\n', rms_error);
+
+ % work out x and y axes
+ ij_corners = [0 0 Nx Nx;
+ 0 Ny Ny 0];
+ ij_tilde = [ij_corners; 1 1 1 1];
+ xy_tilde = P2D^-1 * ij_tilde;
+ lambda = xy_tilde(3,:);
+ xy_corners = [xy_tilde(1,:) ./ lambda;
+ xy_tilde(2,:) ./ lambda];
+ x_lim = sort(xy_corners(1,:), 'ascend');
+ y_lim = sort(xy_corners(2,:), 'ascend');
+ x_lim = x_lim([1 4]);
+ y_lim = y_lim([1 4]);
+
+ x_axis = linspace(x_lim(1), x_lim(2), Nx);
+ y_axis = linspace(y_lim(1), y_lim(2), Ny);
+ i_axis = 1 : Nx;
+ j_axis = 1 : Ny;
+ % calculate dewarping
+ [xmat, ymat] = ndgrid(x_axis, y_axis);
+ ij_tilde_dewarped = P2D * [xmat(:) ymat(:) ones(Nx*Ny,1)]';
+ lambda = ij_tilde_dewarped(3,:);
+ imat = reshape(ij_tilde_dewarped(1,:) ./ lambda, [Nx Ny]);
+ jmat = reshape(ij_tilde_dewarped(2,:) ./ lambda, [Nx Ny]);
+ im_xy = interp2(i_axis, j_axis, im', imat, jmat, 'linear');
+
+ %% show dewarped image
+ figure(1);
+ hold off;
+ imagesc(x_axis, y_axis, im_xy');
+ xlim(x_lim);
+ ylim(y_lim);
+ colormap gray;
+ axis xy equal tight;
+ hold on;
+ set(gca,'Clim', image_clim);
+ for x = (round(x_lim(1)/grid_dx) : round(x_lim(end)/grid_dx)) * grid_dx
+ plot([x x], [y_axis(1) y_axis(end)], 'r-');
+ end
+ for y = (round(y_lim(1)/grid_dy) : round(y_lim(end)/grid_dy)) * grid_dy
+ plot([x_axis(1) x_axis(end)], [y y], 'r-');
+ end
+
+ %% query user to save calibration or repeat
+
+ sInput = input('Would you like to repeat the calibration? y to repeat: ', 's');
+ if strcmpi(sInput, 'y')
+ continue; % continue without incrementing the counter
+ end
+
+ %% save output
+ fprintf('Saving ... ');
+ fstruct = struct('source_image', im, ...
+ 'dewarped_image', im_xy, ...
+ 'plate_pos', plate_origin(:, calib_idx), ...
+ 'P', P2D, ...
+ 'i_axis', i_axis, 'j_axis', j_axis, ...
+ 'x_axis', x_axis, 'y_axis', y_axis, 'x_lim', x_lim, 'y_lim', y_lim, ...
+ 'grid_dx', grid_dx, 'grid_dy', grid_dy, ...
+ 'ij_markers', ij_markers, 'xy_markers', xy_markers, ...
+ 'camera_index', camera_index, ...
+ 'id', camera_id);
+ save(output_file, '-struct', 'fstruct');
+ fprintf('done\n');
+ calib_idx = calib_idx + 1;
+end
\ No newline at end of file
diff --git a/calibration/find_markers.m b/calibration/find_markers.m
new file mode 100644
index 0000000000000000000000000000000000000000..2ed367167736ad29c58151ec9829e82e945f8e17
--- /dev/null
+++ b/calibration/find_markers.m
@@ -0,0 +1,166 @@
+function [ij_markers, xy_markers] = find_markers(img, ij_initial_markers, dot_half_size)
+ % find_markers(img, ij_initial_markers)
+ % conduct a marker search after the user has identified six markers on
+ % the perimeter of a square, arranged like so:
+ %
+ % 5 x x x x 4
+ % x x x x x x
+ % x x x x x x
+ % x x x x x x
+ % 6 x x x x x
+ % 1 2 x x x 3
+
+ %% stuff
+ Nx = size(img, 1);
+ Ny = size(img, 2);
+ thresh = 0.7;
+ grid_tol = 0.1;
+ n_iter = 3;
+
+ %% find x&y position of six markers identified
+ xy_initial_markers= [0 1 nan nan 0 0;
+ 0 0 0 nan nan 1];
+ e_x = ij_initial_markers(:, 3) - ij_initial_markers(:, 1);
+ e_y = ij_initial_markers(:, 5) - ij_initial_markers(:, 1);
+ e_x_norm = e_x / norm(e_x);
+ e_y_norm = e_y / norm(e_y);
+ dx_vec = norm(ij_initial_markers(:, 2) - ij_initial_markers(:,1));
+ dy_vec = norm(ij_initial_markers(:, 6) - ij_initial_markers(:,1));
+
+ xy_initial_markers(1,3) = round( dot(ij_initial_markers(:,3) - ij_initial_markers(:,1), e_x_norm) / dx_vec );
+ xy_initial_markers(2,5) = round( dot(ij_initial_markers(:,5) - ij_initial_markers(:,1), e_y_norm) / dy_vec );
+ xy_initial_markers(1,4) = xy_initial_markers(1,3);
+ xy_initial_markers(2,4) = xy_initial_markers(2,5);
+
+ %% produce initial pinhole model
+ % this maps ij ~ P * xy
+ ij_third = ij_initial_markers(:, 3);
+ ij_initial_markers= ij_initial_markers(:, [1 3 4 5]);
+ xy_initial_markers= xy_initial_markers(:, [1 3 4 5]);
+ P2D = dlt33([ij_initial_markers;ones(1,4)], [xy_initial_markers;ones(1,4)]);
+ for it = 1 : n_iter
+ %% dewarp with initial pinhole model
+ ij_corners = [0 0 Nx Nx;
+ 0 Ny Ny 0];
+ ij_tilde = [ij_corners; 1 1 1 1];
+ xy_tilde = P2D^-1 * ij_tilde;
+ lambda = xy_tilde(3,:);
+ xy_corners = [xy_tilde(1,:) ./ lambda;
+ xy_tilde(2,:) ./ lambda];
+ x_lim = sort(xy_corners(1,:), 'ascend');
+ y_lim = sort(xy_corners(2,:), 'ascend');
+ x_lim = x_lim([1 end]);
+ y_lim = y_lim([1 end]);
+
+ x_axis = linspace(x_lim(1), x_lim(2), Nx);
+ y_axis = linspace(y_lim(1), y_lim(2), Ny);
+ dx = x_axis(2) - x_axis(1);
+ dy = y_axis(2) - y_axis(1);
+
+ [xmat, ymat] = ndgrid(x_axis, y_axis);
+ ij_tilde_dewarped = P2D * [xmat(:) ymat(:) ones(Nx*Ny,1)]';
+ lambda = ij_tilde_dewarped(3,:);
+ i_dewarped = reshape(ij_tilde_dewarped(1,:) ./ lambda, [Nx Ny]);
+ j_dewarped = reshape(ij_tilde_dewarped(2,:) ./ lambda, [Nx Ny]);
+
+ img_dewarped = interp2(1 : Nx, 1 : Ny, img', i_dewarped, j_dewarped);
+
+ %% get dot img
+ [~,ij_third] = min(abs(x_axis));
+ [~,ij_third(2)] = min(abs(y_axis));
+
+ dot_size = dot_half_size * 2 + 1;
+ dot_img_i = round(ij_third(1)) + (-dot_half_size(1) : dot_half_size(1));
+ dot_img_j = round(ij_third(2)) + (-dot_half_size(2) : dot_half_size(2));
+ dot_img = img_dewarped(dot_img_i, dot_img_j);
+
+ %% convolve image and find indices of maxima
+ dot_img = dot_img - mean(dot_img(:));
+ % subtract sliding average so DC component of x-correlation windows is
+ % zero
+ img_sliding_av = filter2(ones(size(dot_img)), img_dewarped) / numel(dot_img);
+ img_sub = img_dewarped - img_sliding_av;
+ % cross-correlate
+ convolved = filter2(dot_img, img_sub);
+ % normalise by image energy
+ image_energy = filter2(ones(size(dot_img)), img_sub.^2);
+ dot_energy = sum(dot_img(:).^2);
+ correl = convolved ./ (sqrt(image_energy) * sqrt(dot_energy));
+ %correl = correl / correl(round(ij_center(1)), round(ij_center(2)));
+
+ correl_left = inf([Nx Ny]);
+ correl_right = inf([Nx Ny]);
+ correl_above = inf([Nx Ny]);
+ correl_below = inf([Nx Ny]);
+
+ correl_left(1:Nx-2, 2:Ny-1) = correl(2:Nx-1, 2:Ny-1);
+ correl_right(3:Nx, 2:Ny-1) = correl(2:Nx-1, 2:Ny-1);
+ correl_above(2:Nx-1, 1:Ny-2) = correl(2:Nx-1, 2:Ny-1);
+ correl_below(2:Nx-1, 3:Ny) = correl(2:Nx-1, 2:Ny-1);
+
+ bw_markers = correl > correl_left ...
+ & correl > correl_right ...
+ & correl > correl_below ...
+ & correl > correl_above ...
+ & correl > thresh;
+
+ %% find xy coordinates of markers
+ [i_markers,j_markers]= find(bw_markers);
+ xy_markers = [x_axis(i_markers(:)'); y_axis(j_markers(:)')];
+ n_markers = size(xy_markers, 2);
+
+ for mark = 1 : n_markers
+ x_range = abs(x_axis - xy_markers(1,mark)) < dx * dot_half_size(1);
+ y_range = abs(y_axis - xy_markers(2,mark)) < dy * dot_half_size(2);
+ sz = [sum(x_range), sum(y_range)];
+ dot_img = img_dewarped(x_range,y_range);
+ dot_img_weighted = dot_img .* fspecial('gaussian', sz, norm(sz)/4);
+ dot_img_x = sum(dot_img_weighted, 2);
+ dot_img_y = sum(dot_img_weighted, 1);
+ xy_markers(1,mark)= sum(dot_img_x' .* x_axis(x_range)) / sum(dot_img_x);
+ xy_markers(2,mark)= sum(dot_img_y .* y_axis(y_range)) / sum(dot_img_y);
+ end
+
+ %% only select markers which are approximately on regular grid
+ xy_frac = xy_markers - round(xy_markers);
+ marker_ok = abs(xy_frac(1,:)) < grid_tol & abs(xy_frac(2,:)) < grid_tol;
+ xy_markers = xy_markers(:, marker_ok);
+ xy_int = round(xy_markers);
+
+ % sometimes, more than one identified marker can be attributed to the same grid
+ % point. eliminate these, based on which marker is closest
+ [~,ind] = unique(xy_int', 'rows');
+ n_unique = length(ind);
+ xy_markers_uniq = zeros(2, n_unique);
+ for n = 1 : n_unique
+ ndist = sum(bsxfun(@minus, xy_markers,xy_int(:,ind(n))).^2, 1);
+ [~,inearest] = min(ndist);
+ xy_markers_uniq(:,n) = xy_markers(:, inearest);
+ end
+ xy_markers = xy_markers_uniq;
+ n_markers = n_unique;
+
+ %% get i,j coord of markers
+ xy_tilde = [xy_markers; ones(1, n_markers)];
+ ij_tilde = P2D * xy_tilde;
+ lambda = [1;1] * ij_tilde(3,:);
+ ij_markers = ij_tilde(1:2,:) ./ lambda;
+
+ %% calculate a new calibration
+ xy_markers = round(xy_markers);
+ P2D = dlt33([ij_markers;ones(1,n_markers)], [xy_markers; ones(1,n_markers)]);
+ end
+
+ %% display
+ figure(10);
+ hold off;
+ imagesc(1 : Nx, 1 : Ny, img');
+ xlim([1 Nx]);
+ ylim([1 Ny]);
+ axis equal tight xy;
+ colormap gray;
+ hold on;
+ plot(ij_markers(1,:), ij_markers(2,:), 'r.');
+end
+
+
diff --git a/calibration/simple_calibration.m b/calibration/simple_calibration.m
new file mode 100644
index 0000000000000000000000000000000000000000..7e97a2e99b602e8d0fea19529036e2a62f8187db
--- /dev/null
+++ b/calibration/simple_calibration.m
@@ -0,0 +1,223 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% bulk_marker_search.m
+% find markers in images of calibration plates
+% then perform a 2D calibration to dewarp the calibration plate image
+
+addpath([pwd '/../../mex/minmax/']);
+addpath([pwd '/dlt/']);
+addpath([pwd '/poly/']);
+addpath([pwd '/../pinhole/']);
+%% input parameters
+
+grid_dx = 20; % real units, mm
+grid_dy = 20;
+invert_image = 1; % invert image for black dots on white
+blur_image = 0; % gaussian blur image
+LoG_image = 0; % apply laplacian of Gaussian filter
+normalise_image = 0;
+
+
+output_dir = 'M:\mixer\calibration\demo\';
+input_file_list = { 'M:\mixer\calibration\demo\B00001.tif'};
+plate_origin = [0 0 0]';
+camera_index = 1;
+camera_id = 1;
+dot_half_window = [32 32];
+calib_idx_offset = 0;
+origin_pos = [0 0];
+
+image_clim = [0 255];
+
+%% generate calibration for each image
+
+nCalibrations = length(input_file_list);
+
+calib_idx = 1;
+
+while calib_idx <= nCalibrations
+ %% update user of status
+ input_file = input_file_list{calib_idx};
+ output_file = sprintf('%s/CAM%d-%02d.mat', output_dir, camera_index, calib_idx+calib_idx_offset);
+ fprintf('Image %02d of %02d:\n', calib_idx, nCalibrations);
+
+ % check if calibration already exists
+ if exist(output_file, 'file')
+ sInput = input('Calibration already exists - overwrite? y to overwrite: ', 's');
+ if ~strcmpi(sInput, 'y')
+ calib_idx = calib_idx + 1;
+ continue;
+ end
+ end
+ % check if image input exists
+ if ~exist(input_file, 'file')
+ fprintf(' - could not find image file\n');
+ end
+
+ %% load image, request position
+ im = imread(input_file, 2)';
+ im = double(im);
+ im = fliplr(im);
+ Nx = size(im, 1);
+ Ny = size(im, 2);
+
+ %% create filtered image for processing
+ im_filt = im;
+ if invert_image
+ im_filt = max(im(:))-im;
+ end
+ if blur_image
+ G = fspecial('gaussian', [5 5]);
+ im_filt = filter2(G, im_filt);
+ end
+ if LoG_image
+ scale = norm(dot_half_window)*0.2;
+ G = fspecial('log', 2*dot_half_window, scale);
+ im_filt = scale * conv2(G, im_filt);
+ end
+ if normalise_image
+ [minv, maxv] = minmaxfiltnd(single(im_filt), dot_half_window*2+1);
+ im_filt = (im_filt - minv) ./ max(16, maxv-minv);
+ end
+
+ %% display image
+ figure(1);
+ hold off;
+ imagesc(im_filt');
+ axis equal tight xy;
+ colormap gray;
+ %set(gca,'Clim', image_clim);
+ hold on;
+ colorbar;
+
+ %% select region of interest
+ fprintf('Please select region of interest');
+ [i_coord, j_coord] = ginput(2);
+ xlim(i_coord);
+ ylim(j_coord);
+
+ %% get six marker positions
+ % center, first right, far right, far corner, far above, above
+ bDone = false;
+ while ~bDone
+ %% get marker position from user
+ fprintf('Please identify the six markers\n');
+ fprintf('Press any key to abort and try again\n');
+
+ i_marker = nan(1,6);
+ j_marker = nan(1,6);
+ b_aborted = false;
+ h_plot = [];
+ for mark = 1 : 6
+ %% get marker position from user input
+ [i_marker(mark) j_marker(mark), btn] = ginput(1);
+ if btn > 3
+ b_aborted= true;
+ break;
+ end
+ %% get barycentre of marker, to be more accurate
+ dot_img_i = round(i_marker(mark)) + (-dot_half_window(1) : dot_half_window(1));
+ dot_img_j = round(j_marker(mark)) + (-dot_half_window(2) : dot_half_window(2));
+ [imat,jmat] = ndgrid(dot_img_i, dot_img_j);
+ dot_img = im_filt(dot_img_i, dot_img_j);
+ i_marker(mark) = sum(dot_img(:) .* imat(:)) / sum(dot_img(:));
+ j_marker(mark) = sum(dot_img(:) .* jmat(:)) / sum(dot_img(:));
+
+ %% show user what he's plotting
+ delete(h_plot);
+ h_plot = plot(i_marker, j_marker, 'ro-');
+ for it = 1 : mark
+ plot( i_marker(it)+[-1 -1 +1 +1 -1]*dot_half_window(1), ...
+ j_marker(it)+[-1 +1 +1 -1 -1]*dot_half_window(2), 'g-');
+ end
+ end
+ if b_aborted
+ delete(h_plot);
+ continue;
+ end
+ bDone = true;
+ end
+
+ %% find marker positions
+ [ij_markers, xy_markers]= find_markers(im_filt, [i_marker;j_marker], dot_half_window);
+ xy_markers(1,:) = xy_markers(1,:) * grid_dx + origin_pos(1);
+ xy_markers(2,:) = xy_markers(2,:) * grid_dy + origin_pos(2);
+ n_markers = size(xy_markers,2);
+ xy_tilde = [xy_markers; ones(1,n_markers)];
+
+ %% create calibration and dewarp
+ P2D = dlt33([ij_markers;ones(1,n_markers)], xy_tilde);
+ ij_tilde = P2D * xy_tilde;
+ lambda = [1; 1] * ij_tilde(3,:);
+ ij_fit = ij_tilde(1:2,:) ./ lambda;
+
+ ij_error = ij_markers - ij_fit;
+ px_error = sqrt(ij_error(1,:).^2 + ij_error(2,:).^2);
+ rms_error = sqrt(mean(px_error.^2));
+
+ fprintf('RMS fit error: %0.2f px\n', rms_error);
+
+ % work out x and y axes
+ ij_corners = [0 0 Nx Nx;
+ 0 Ny Ny 0];
+ ij_tilde = [ij_corners; 1 1 1 1];
+ xy_tilde = P2D^-1 * ij_tilde;
+ lambda = xy_tilde(3,:);
+ xy_corners = [xy_tilde(1,:) ./ lambda;
+ xy_tilde(2,:) ./ lambda];
+ x_lim = sort(xy_corners(1,:), 'ascend');
+ y_lim = sort(xy_corners(2,:), 'ascend');
+ x_lim = x_lim([1 4]);
+ y_lim = y_lim([1 4]);
+
+ x_axis = linspace(x_lim(1), x_lim(2), Nx);
+ y_axis = linspace(y_lim(1), y_lim(2), Ny);
+ i_axis = 1 : Nx;
+ j_axis = 1 : Ny;
+ % calculate dewarping
+ [xmat, ymat] = ndgrid(x_axis, y_axis);
+ ij_tilde_dewarped = P2D * [xmat(:) ymat(:) ones(Nx*Ny,1)]';
+ lambda = ij_tilde_dewarped(3,:);
+ imat = reshape(ij_tilde_dewarped(1,:) ./ lambda, [Nx Ny]);
+ jmat = reshape(ij_tilde_dewarped(2,:) ./ lambda, [Nx Ny]);
+ im_xy = interp2(i_axis, j_axis, im', imat, jmat, 'linear');
+
+ %% show dewarped image
+ figure(1);
+ hold off;
+ imagesc(x_axis, y_axis, im_xy');
+ xlim(x_lim);
+ ylim(y_lim);
+ colormap gray;
+ axis xy equal tight;
+ hold on;
+ set(gca,'Clim', image_clim);
+ for x = (round(x_lim(1)/grid_dx) : round(x_lim(end)/grid_dx)) * grid_dx
+ plot([x x], [y_axis(1) y_axis(end)], 'r-');
+ end
+ for y = (round(y_lim(1)/grid_dy) : round(y_lim(end)/grid_dy)) * grid_dy
+ plot([x_axis(1) x_axis(end)], [y y], 'r-');
+ end
+
+ %% query user to save calibration or repeat
+
+ sInput = input('Would you like to repeat the calibration? y to repeat: ', 's');
+ if strcmpi(sInput, 'y')
+ continue; % continue without incrementing the counter
+ end
+
+ %% save output
+ fprintf('Saving ... ');
+ fstruct = struct('source_image', im, ...
+ 'dewarped_image', im_xy, ...
+ 'plate_pos', plate_origin(:, calib_idx), ...
+ 'P', P2D, ...
+ 'i_axis', i_axis, 'j_axis', j_axis, ...
+ 'x_axis', x_axis, 'y_axis', y_axis, 'x_lim', x_lim, 'y_lim', y_lim, ...
+ 'grid_dx', grid_dx, 'grid_dy', grid_dy, ...
+ 'ij_markers', ij_markers, 'xy_markers', xy_markers, ...
+ 'camera_index', camera_index, ...
+ 'id', camera_id);
+ save(output_file, '-struct', 'fstruct');
+ fprintf('done\n');
+ calib_idx = calib_idx + 1;
+end
\ No newline at end of file
diff --git a/image_preprocess.m b/image_preprocess.m
index 42fe1a2650580a88f161f74469abb0c4c4ea0abe..e4336ccf3a5e694a3b66b0cb76ef71f47dc57a31 100644
--- a/image_preprocess.m
+++ b/image_preprocess.m
@@ -37,20 +37,46 @@ function im_filt = image_preprocess(im, filters)
if isfield(filter, 'size'), sz = filter.size; end
if isfield(filter, 'sigma'), sigma = filter.sigma; end
G = fspecial('gaussian', sz, sigma);
- im_filt = filter2(G, im_filt);
+ im_filt = filter2(G, im_filt);
+ elseif strcmp(filter.type, 'custom')
+ %% custom filter specified with fspecial
+ G = fspecial('gaussian', [3 3]);
+ if isfield(filter, 'kernel'), G = filter.kernel; end
+ im_filt = filter2(G, im_filt);
elseif strcmp(filter.type, 'ssmin')
%% subtract sliding minimum
sz = [7 7];
if isfield(filter, 'size'), sz = filter.size; end
sz = sz + mod(sz+1,2);
-
- minval = minmaxfiltnd(im_filt, sz);
- im_filt = im_filt - minval;
+ minval = medfilt2(im_filt, [3 3]); % median filter first to avoid noise
+ % sliding minimum filter to get local minimum value
+ minval = minfiltnd(minval, sz);
+ % box filter to smooth background image
+ minval = convn(minval, ones(1, sz(1))/sz(1), 'same');
+ minval = convn(minval, ones(sz(2), 1)/sz(2), 'same');
+ % clipping operation to enforce positivity
+ im_filt = max(0, im_filt - minval);
+ elseif strcmp(filter.type, 'lmax')
+ %% local maximum
+ sz = [7 7];
+ if isfield(filter, 'size'), sz = filter.size; end
+ sz = sz + mod(sz+1,2);
+ im_filt = maxfiltnd(im_filt, sz);
elseif strcmp(filter.type, 'ssbg')
%% subtract background
- bg = zeros(size(im));
+ bg = zeros(size(im));
if isfield(filter, 'bg'), bg = filter.bg; end
- im_filt = max(0, im - bg);
+ im_filt = max(0, im - bg);
+ elseif strcmp(filter.type, 'levelize')
+ %% levelize
+ white = ones(size(im));
+ if isfield(filter, 'white'), white = filter.white; end
+ im_filt = im_filt ./ white;
+ elseif strcmp(filter.type, 'median')
+ %% median filtering
+ sz = [5 5];
+ if isfield(filter, 'size'), sz = filter.size; end
+ im_filt = medfilt2(im_filt, sz, 'symmetric');
elseif strcmp(filter.type, 'norm')
%% normalisation filter
sz = [7 7];
@@ -61,15 +87,60 @@ function im_filt = image_preprocess(im, filters)
sz = sz + mod(sz+1,2);
[minv, maxv] = minmaxfiltnd(im_filt, sz);
- im_filt = (im_filt - minv) ./ max(maxv-minv, 1/max_gain);
+ im_filt = (im_filt - minv) ./ max(maxv-minv, 1/max_gain);
+ elseif strcmp(filter.type, 'norm2')
+ %% normalisation filter
+ sz = [7 7];
+ max_gain = 1;
+
+ if isfield(filter, 'size'), sz = filter.size; end
+ if isfield(filter, 'max_gain'), max_gain = filter.max_gain; end
+ sz = sz + mod(sz+1,2);
+
+ [minv, maxv] = minmaxfiltnd(im_filt, sz);
+ % box filter local maximum
+ maxv = convn(maxv, ones(1, sz(1))/sz(1), 'same');
+ maxv = convn(maxv, ones(sz(2), 1)/sz(2), 'same');
+ % box filter local minimum
+ minv = convn(minv, ones(1, sz(1))/sz(1), 'same');
+ minv = convn(minv, ones(sz(2), 1)/sz(2), 'same');
+ im_filt = (im_filt - minv) ./ max(maxv-minv, 1/max_gain);
+ elseif strcmp(filter.type, 'mednorm')
+ %% normalisation filter with median filtering
+ sz = [7 7];
+ max_gain = 1;
+
+ if isfield(filter, 'size'), sz = filter.size; end
+ if isfield(filter, 'max_gain'), max_gain = filter.max_gain; end
+ sz = sz + mod(sz+1,2);
+
+ % local maximum
+ maxv = minmaxfiltnd(im_filt, sz);
+ % box filter local maximum
+ maxv = convn(maxv, ones(1, sz(1))/sz(1), 'same');
+ maxv = convn(maxv, ones(sz(2), 1)/sz(2), 'same');
+ im_filt = max(0,im_filt) ./ max(maxv, 1/max_gain);
elseif strcmp(filter.type, 'clip')
- %% clip filter
- im_med = median(im_filt(:));
- im_std = std(im_filt(:));
- n = 2;
- if isfield(filter, 'n'), n = filter.n; end
- thresh = im_med + n*im_std;
- im_filt = min(im_filt, thresh);
+ %% clip filter, with automatic
+ if isfield(filter, 'threshold')
+ % manually specified threshold
+ threshold = filter.threshold;
+ im_filt = max(threshold(1), min(threshold(2), im_filt));
+ else
+ % algorithmically specified threshold
+ im_med = median(im_filt(:));
+ im_std = std(im_filt(:));
+ n = 2;
+ if isfield(filter, 'n'), n = filter.n; end
+ thresh = im_med + n*im_std;
+ im_filt = min(im_filt, thresh);
+ end
+ elseif strcmp(filter.type, 'transpose')
+ im_filt = im_filt.';
+ elseif strcmp(filter.type, 'invert')
+ offset = 0;
+ if isfield(filter, 'offset'), offset = filter.offset; end
+ im_filt = offset - im_filt;
elseif strcmp(filter.type, 'null')
%% null filter, do nothing
im_filt = im_filt;
diff --git a/job_create.m b/job_create.m
index 199fccb6301079d25c2e263593be21954640e0a6..d8b11bd6fea492e4f210c35cd5c2e7709c0e3be0 100644
--- a/job_create.m
+++ b/job_create.m
@@ -10,7 +10,7 @@ addpath(genpath(fileparts(mfilename('fullpath'))));
%% input parameters: file paths
% where to save jobs
job_idx = 1;
-job_dir = 'J:\FEE\EFMRL\John\TR-IM7\jobs\';
+job_dir = 'F:\mixer\0609\run01\';
job_pat = '%s/job%03d.mat';
% timeseries PIV flag
@@ -18,16 +18,16 @@ job_pat = '%s/job%03d.mat';
% in a timeseries
b_timeseries = true;
% path of images to process
-src_dir = 'J:\FEE\EFMRL\John\TR-IM7';
-src_pat = '%s/B%05d.im7'; % source image to process
-src_range = 1:3; % range of source images to process
+src_dir = 'M:\mixer\151020\dt1.00\';
+src_pat = '%s/B%05d.tif'; % source image to process
+src_range = 1:10; % range of source images to process
% background & mask
-bg_dir = src_dir; % background image directory
-bg_pat = '%s/minmax.mat'; % background image (ignore if not using)
+bg_dir = ''; % background image directory
+bg_pat = ''; % background image (ignore if not using)
msk_dir = [src_dir '/MASK/']; % mask image directory
msk_pat = '%s/B0002.im7'; % pattern for mask image input (ignore if not using)
% where to save result
-vec_dir = [src_dir '/C12']; % directory to save vector fields to
+vec_dir = [src_dir '/C1']; % directory to save vector fields to
vec_pat = '%s/B%05d.mat'; % output vector field file pattern
%% input parameters: processing
@@ -37,7 +37,7 @@ b_overwrite = true; % overwrite existing vector fields?
b_execute_now = true; % start processing now?
% PIV correlation settings
-i_frames = {[1 1],[2 2], [3 3]}; % indices of frame pairs to cross-correlate
+i_frames = {[1 2]}; % indices of frame pairs to cross-correlate
n_passes = 3; % number of passes
n_peaks = 3; % number of correlation peaks to detect
wsize = [64 64; 32 32; 32 32];% window size
@@ -122,6 +122,7 @@ for n = 1 : n_src
end
end
fprintf('Loading sample %s\n', src_name);
+
% load image
imx = readimx(src_name);
im_ary = cellfun(@(f) double(f.Components{1}.Planes{1}), imx.Frames, 'UniformOutput', false);
diff --git a/job_create_2d.m b/job_create_2d.m
new file mode 100644
index 0000000000000000000000000000000000000000..abd6474e3de42f2d24e54c525fda80ff722da877
--- /dev/null
+++ b/job_create_2d.m
@@ -0,0 +1,270 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% create job file for PIV processing
+%
+% Author: John Lawson
+% j.m.lawson@soton.ac.uk
+% Created: 11/09/2019
+%
+addpath(genpath(fileparts(mfilename('fullpath'))));
+addpath([pwd '/../misc']);
+
+%% input parameters: file paths
+% where to save jobs
+job_idx = 24;
+job_dir = 'sampledata';
+job_pat = '%s/job%03d.mat';
+
+% timeseries PIV flag
+% in timeseries PIV, correlations are conducted between subsequent images
+% in a timeseries
+b_timeseries = true;
+% path of images to process
+dset = sprintf('run%02d', job_idx);
+src_dir = ['sampledata/'];
+src_pat = '%s/NACA0020_%05d.jpg'; % source image to process
+src_range = 0:1; % range of source images to process
+% background & mask
+bg_dir = src_dir; % background image directory
+bg_pat = '%s/Bmin.tif'; % background image (ignore if not using)
+msk_dir = [src_dir '/MASK/']; % mask image directory
+msk_pat = '%s/B0002.im7'; % pattern for mask image input (ignore if not using)
+% where to save result
+vec_dir = ['sampledata/']; % directory to save vector fields to
+vec_pat = '%s/B%05d.mat'; % output vector field file pattern
+% calibration
+b_transpose = true;
+b_apply_calib = false;
+cal_dir = 'M:\mixer\calibration\20-12-09_02\';
+cal_pat = '%s/CAM1-01.mat';
+
+%% input parameters: processing
+% processing flags
+b_use_mask = false; % use mask file?
+b_overwrite = true; % overwrite existing vector fields?
+b_execute_now = true; % start processing now?
+
+% PIV correlation settings
+i_frames = {[1 1]}; % indices of frame pairs to cross-correlate
+n_passes = 3; % number of passes
+n_peaks = 3; % number of correlation peaks to detect
+wsize = [64 64; 32 32; 32 32];% window size
+wtype = {'gaussian','gaussian','gaussian','gaussian'}; % window weighting function type
+overlap = [50 50 50]; % window overlap %
+peak_finder = 'gauss3'; % peak finder:
+ % gauss 3 (3 point fit)
+ % gauss 4 (4 parameter lsq gaussian fit)
+ % gauss 5 (5 parameter lsq gaussian fit)
+ % gauss 6 (6 parameter lsq gaussian fit)
+
+%% input parameters: plotting features
+plt_opts = struct();
+plt_opts.c_lim = [0 1]*256;
+plt_opts.font_size = 24; % figure font size
+plt_opts.b_plot_src = true & b_execute_now; % plot source images?
+plt_opts.b_plot_u = false & b_execute_now; % plot u component field?
+plt_opts.b_plot_v = false & b_execute_now; % plot v component field?
+plt_opts.b_plot_cv = false & b_execute_now; % plot correlation value?
+
+%% input parameters: image preprocessing toolchain
+% add filters in order in which they are to be applied
+% see image_preprocess for more details and options of filters
+min_filt_size = [7 7];% [11 11]; % sliding bg filter size
+gauss_filt_size = [3 3]; % gaussian filter size
+norm_filt_size = [9 9]; % normalisation filter size
+log_filt_size = [9 9]; % Laplacian-of-gaussian filter size
+log_sigma = 1.2; % Laplacian-of-Gaussian filter width
+
+filters = {struct('type', 'null')};
+
+% % subtract background image: will require background image to be supplied
+%filters{end+1} = struct('type', 'ssbg');
+% % subtract sliding minimum
+%filters{end+1} = struct('type', 'ssmin', 'size', min_filt_size);
+% % max-clip bright particles at threshold = median(im) + n*sigma(im)
+filters{end+1} = struct('type', 'clip', 'n', 2);
+% % 2D Laplacian of Gaussian filter
+%G_filt = -fspecial('log', log_filt_size, log_sigma);
+%filters{end+1} = struct('type', 'custom', 'kernel', G_filt);
+% 2D gaussian filter
+filters{end+1} = struct('type', 'gaussian', 'size', gauss_filt_size, 'sigma', 0.67);
+% % min-max normalisation filter
+%filters{end+1} = struct('type', 'norm2', 'size', norm_filt_size, 'max_gain', 1/32);
+%filters{end+1} = struct('type', 'mednorm', 'size', norm_filt_size, 'max_gain', 1);
+% transpose image
+if b_transpose
+ filters{end+1} = struct('type', 'transpose');
+end
+
+%% input parameters: vector postprocesing
+% load calibration
+operations = {struct('type', 'null')};
+
+if b_apply_calib
+ cal_file = sprintf(cal_pat,cal_dir);
+ calib = matfile(cal_file,'writable',false);
+ P2D = calib.P;
+
+ operations{end+1} = struct('type', 'scale2d', ...
+ 'P', P2D);
+end
+
+%% check job directory
+job_completed_dir = [job_dir '/completed/'];
+job_busy_dir = [job_dir '/busy/'];
+job_file = sprintf(job_pat, job_dir, job_idx);
+
+fprintf('------------------------------------------------------------------\n');
+fprintf('Creating job %s\n', job_file);
+% check directory exists
+if ~exist(job_dir, 'dir')
+ mkdir(job_dir);
+ %error('job directory %s does not exist');
+end
+% check complete and busy directories are there
+if ~exist(job_completed_dir,'dir')
+ fprintf('making completed jobs directory ...\n');
+ mkdir(job_completed_dir);
+end
+if ~exist(job_busy_dir,'dir')
+ fprintf('making busy jobs directory ...\n');
+ mkdir(job_busy_dir);
+end
+% check if existing job file already exists
+if exist(job_file, 'file')
+ s_override = input('job already exists, are you sure you want to overwrite? y/N ', 's');
+ if ~strcmp(s_override, 'y')
+ error('will not overwrite job');
+ end
+end
+% make output directory
+if ~exist(vec_dir, 'dir')
+ fprintf('making output directory ...\n');
+ mkdir(vec_dir);
+end
+
+%% load sample and check all files in range exist
+n_src = length(src_range);
+fprintf('------------------------------------------------------------------\n');
+fprintf('Checking %d files in %s ...\n', n_src, src_dir);
+for n = 1 : n_src
+ src_name = sprintf(src_pat, src_dir, src_range(n));
+ if ~exist(src_name, 'file')
+ error('source %s does not exist\n', src_name);
+ end
+end
+fprintf('Loading sample %s\n', src_name);
+% load image
+im_ary = readimg(src_name);
+% image properties
+n_frames = length(im_ary);
+im_sz = size(im_ary{1});
+assert(all(all(cell2mat(i_frames)>=1 & cell2mat(i_frames) <= n_frames)), 'PIV correlations must be specified between valid frame indices');
+
+% report
+fprintf('frames: %d\n', n_frames);
+fprintf('resolution: %d x %d\n', im_sz);
+
+%% check mask
+if b_use_mask
+ msk_name = sprintf(msk_pat, msk_dir);
+ fprintf('Checking mask %s ...\n', msk_name);
+ im_mask_ary = readimg(msk_name);
+ assert(length(im_mask_ary) == n_frames, 'mask image does not contain the right number of frames');
+ assert(all(size(im_mask_ary{1}) == im_sz), 'mask image is not the same size as source');
+end
+
+%% load background image
+% flag to indicate we need to load
+b_use_bg = any(cellfun(@(f)strcmp(f.type,'ssbg'),filters));
+
+if b_use_bg
+ bg_name = sprintf(bg_pat, bg_dir);
+ fprintf('Checking background %s\n', bg_name);
+ assert(~~exist(bg_name,'file'), 'background image is required, but file is missing');
+ im_bg_ary = readimg(bg_name);
+ assert(length(im_bg_ary) == n_frames, 'background image does not contain the right number of frames');
+ assert(all(size(im_bg_ary{1}) == im_sz), 'background image is not the same size as source');
+end
+
+%% define PIV image grid
+% zero indexed
+% TODO: make this match the internal scale/offset of im7 image files
+if b_transpose
+ Nx = im_sz(2);
+ Ny = im_sz(1);
+else
+ Nx = im_sz(1);
+ Ny = im_sz(2);
+end
+im_x = 0 : Nx-1;
+im_y = 0 : Ny-1;
+
+%% define PIV options structure
+% create piv options array
+% note clever trick with i_frames cell array, which will generate array of
+% options
+piv_opts = struct('im_x', im_x, ...
+ 'im_y', im_y, ...
+ 'dt', 1, ...
+ 'i_frames', i_frames, ...
+ 'n_passes', n_passes, ...
+ 'n_peaks', n_peaks, ...
+ 'wsize', wsize, ...
+ 'overlap', overlap, ...
+ 'peak_finder', peak_finder, ...
+ 'median_filter', true, ...
+ 'epsilon_0', 0.2, ... % typical PIV error, px
+ 'epsilon_thresh', 2, ... % replace if vector exceeds n*RMS of neighbours
+ 'mask_thresh', 0.25); % maximum acceptable fraction of masked pixels in IW
+% add window type field
+n_correl = length(piv_opts);
+for n = 1 : n_correl
+ piv_opts(n).wtype = wtype;
+end
+
+%% bundle into structure
+job = struct();
+% job details
+job.job_idx = job_idx;
+job.job_dir = job_dir;
+job.job_pat = job_pat;
+% timeseries flag
+job.b_timeseries= b_timeseries;
+% source images
+job.src_dir = src_dir;
+job.src_pat = src_pat;
+job.src_range = src_range;
+% background and mask
+job.b_use_bg = b_use_bg;
+job.bg_dir = bg_dir;
+job.bg_pat = bg_pat;
+job.b_use_mask = b_use_mask;
+job.msk_dir = msk_dir;
+job.msk_pat = msk_pat;
+% output vector fields
+job.b_overwrite=b_overwrite;
+job.vec_dir = vec_dir;
+job.vec_pat = vec_pat;
+
+% filtering
+job.filters = filters;
+
+% PIV processing options
+job.n_correl = n_correl;
+job.piv_opts = piv_opts;
+
+% postprocessing operations
+job.operations = operations;
+
+% plotting options
+job.plt_opts = plt_opts;
+
+%% save job to file
+fprintf('saving job ... \n');
+save(job_file, '-struct', 'job');
+
+%% execute
+if b_execute_now
+ job_execute(job_file);
+ delete(job_file);
+end
\ No newline at end of file
diff --git a/job_create_postproc.m b/job_create_postproc.m
new file mode 100644
index 0000000000000000000000000000000000000000..3dd5cbe074f964be3025b0d7d056be89d6a13b40
--- /dev/null
+++ b/job_create_postproc.m
@@ -0,0 +1,108 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% create job file for PIV processing
+%
+% Author: John Lawson
+% j.m.lawson@soton.ac.uk
+% Created: 11/09/2019
+%
+addpath(genpath(fileparts(mfilename('fullpath'))));
+
+%% input parameters: file paths
+% where to save jobs
+job_idx = 1;
+job_dir = 'F:\mixer\091020\jobs\';
+job_pat = '%s/pproc%03d.mat';
+
+% path of vector fields to postprocess
+src_dir = ['F:\mixer\091020\' sprintf('run%02d', job_idx)];
+src_pat = '%s/C1/B%05d.mat'; % source image to process
+src_range = 1:1000; % range of source images to process
+% where to save result
+dst_dir = src_dir; % directory to save vector fields to
+dst_pat = '%s/V%05d.mat'; % output vector field file pattern
+% 2D calibration
+cal_dir = 'F:\mixer\061020\calibration\';
+cal_pat = '%s/CAM1-01.mat';
+piv_dt_range = [20 20 20 20 15 7.5 30 20 15]*1E-3;
+
+%% input parameters: vector postprocesing
+% load calibration
+cal_file = sprintf(cal_pat,cal_dir);
+calib = load(cal_file);
+P2D = calib.P;
+
+operations = {struct('type', 'null')};
+operations{end+1} = struct('type', 'scaledt', ...
+ 'dt', piv_dt_range(job_idx));
+operations{end+1} = struct('type', 'scale2d', ...
+ 'P', P2D);
+
+%% check job directory
+job_completed_dir = [job_dir '/completed/'];
+job_busy_dir = [job_dir '/busy/'];
+job_file = sprintf(job_pat, job_dir, job_idx);
+
+fprintf('------------------------------------------------------------------\n');
+fprintf('Creating job %s\n', job_file);
+% check directory exists
+if ~exist(job_dir, 'dir')
+ error('job directory %s does not exist');
+end
+% check complete and busy directories are there
+if ~exist(job_completed_dir,'dir')
+ fprintf('making completed jobs directory ...\n');
+ mkdir(job_completed_dir);
+end
+if ~exist(job_busy_dir,'dir')
+ fprintf('making busy jobs directory ...\n');
+ mkdir(job_busy_dir);
+end
+% check if existing job file already exists
+if exist(job_file, 'file')
+ s_override = input('job already exists, are you sure you want to overwrite? y/N ', 's');
+ if ~strcmp(s_override, 'y')
+ error('will not overwrite job');
+ end
+end
+% make output directory
+if ~exist(dst_dir, 'dir')
+ fprintf('making output directory ...\n');
+ mkdir(dst_dir);
+end
+
+%% load sample and check all files in range exist
+n_src = length(src_range);
+fprintf('------------------------------------------------------------------\n');
+fprintf('Checking %d files in %s ...\n', n_src, src_dir);
+for n = 1 : n_src
+ src_name = sprintf(src_pat, src_dir, src_range(n));
+ if ~exist(src_name, 'file')
+ error('source %s does not exist\n', src_name);
+ end
+end
+
+%% bundle into structure
+job = struct();
+% job details
+job.job_idx = job_idx;
+job.job_dir = job_dir;
+job.job_pat = job_pat;
+% source images
+job.src_dir = src_dir;
+job.src_pat = src_pat;
+job.src_range = src_range;
+% output vector fields
+job.b_overwrite=b_overwrite;
+job.dst_dir = dst_dir;
+job.dst_pat = dst_pat;
+% postprocessing operations
+job.operations = operations;
+
+%% save job to file
+fprintf('saving job ... \n');
+save(job_file, '-struct', 'job');
+
+%% execute
+if b_execute_now
+ job_postprocess(job_file);
+end
\ No newline at end of file
diff --git a/job_execute.m b/job_execute.m
index 6d1d3d6f2718fc7f27f2b19d6be11e812561b8af..ffef2fe5fea6379e22353047a0c1c13e5d0c3ca0 100644
--- a/job_execute.m
+++ b/job_execute.m
@@ -30,6 +30,11 @@ function job_execute(job_file)
piv_opts = job.piv_opts;
% plotting options
plt_opts = job.plt_opts;
+ % postprocessing
+ operations = {};
+ if isfield(job, 'operations')
+ operations = job.operations;
+ end
%% report
fprintf('--------------------------------------------------------------\n');
@@ -39,8 +44,10 @@ function job_execute(job_file)
%% load prototype image
% load first image
src_name = sprintf(src_pat, src_dir, src_range(1));
- imx = readimx(src_name);
- im_ary = cellfun(@(f) double(f.Components{1}.Planes{1}), imx.Frames, 'UniformOutput', false);
+ im_ary = readimg(src_name);
+ % apply test preprocessing to first image
+ % this handles things like changing image size
+ im_ary{1} = image_preprocess(im_ary{1}, filters);
% extract properties
n_frames = length(im_ary);
im_sz = size(im_ary{1});
@@ -49,8 +56,9 @@ function job_execute(job_file)
if b_use_bg
bg_name = sprintf(bg_pat, bg_dir);
fprintf('Load background %s\n', bg_name);
- s_bg = importdata(bg_name);
- im_bg_ary = s_bg.min;
+ %s_bg = importdata(bg_name);
+ %im_bg_ary = s_bg.min;
+ im_bg_ary = readimg(bg_name);
else
im_bg_ary = {zeros(im_sz)};
im_bg_ary(1:n_frames) = im_bg_ary;
@@ -60,8 +68,7 @@ function job_execute(job_file)
if b_use_mask
msk_name = sprintf(msk_pat, msk_dir);
fprintf('Load mask %s\n', msk_name);
- imx_mask = readimx(msk_name);
- im_mask_ary = cellfun(@(f) f.Components{2}.Planes{1} == 0, imx_mask.Frames, 'UniformOutput', false);
+ im_mask_ary = readimg(msk_name);
else
im_mask_ary = {false(im_sz)};
im_mask_ary(1:n_frames) = im_mask_ary;
@@ -85,19 +92,19 @@ function job_execute(job_file)
continue;
end
- % load image
- imx = readimx(src_name);
- im_ary = cellfun(@(f) double(f.Components{1}.Planes{1}), imx.Frames, 'UniformOutput', false);
-
- % read attributes
- atts = readAttributes(imx);
- if isKey(atts, 'Reference time dt 1')
- s = atts('Reference time dt 1');
- piv_dt = 1E-6 * s.value;
- else
- warning('dt field is missing from IM7 file, dt = 1');
- piv_dt = 1;
- end
+ % load image and attributes
+ [im_ary, atts] = readimg(src_name);
+
+ %% read attributes
+ if isKey(atts, 'Reference time dt 1')
+ s = atts('Reference time dt 1');
+ piv_dt = 1E-6 * s.value;
+ elseif isKey(atts, 'DeltaT')
+ piv_dt = atts('DeltaT');
+ else
+ warning('dt field is missing from file, dt = 1');
+ piv_dt = 1;
+ end
t_load = toc(t_samp);
%% preprocess images
@@ -155,34 +162,57 @@ function job_execute(job_file)
im_A = im_pairs{1,c};
im_B = im_pairs{2,c};
+
+
+ %% cross-correlate
+ fprintf(' correlation %d of %d\n', c, n_correl);
+ piv_result{c} = PIV_2D_wdef(im_A, im_B, ...
+ im_correl_mask{c}, ...
+ piv_opts(c));
+
%% plot source, if requested
if plt_opts.b_plot_src
+ %% result at final pass
+ V = piv_result{c}(end);
+
+ %% source image
figure(10);
clf;
- imshowpair(im_A', im_B');
+ ax_src = gca;
+ h_pair = imshowpair(im_A', im_B');
+ h_pair.XData = V.win_x;
+ h_pair.YData = V.win_y;
axis equal tight ij;
+
+ %% overlay vector field
+ V.peak_mag(V.nan_mask)= nan;
+ hold on;
+ [xm,ym] = ndgrid(V.win_x, V.win_y);
+ ir = 1 : 2 : V.n_windows(1);
+ jr = 1 : 2 : V.n_windows(2);
+ quiversc(xm(ir,jr),ym(ir,jr),V.ux(ir,jr),V.uy(ir,jr),10,'k');
+
drawnow;
end
-
- %% cross-correlate
- fprintf(' correlation %d of %d\n', c, n_correl);
- piv_result{c} = PIV_2D_wdef(im_A, im_B, ...
- im_correl_mask{c}, ...
- piv_opts(c));
end
t_proc = toc(t_proc);
-
- %% save result
- t_save = tic;
- fprintf(' save to %s\n', vec_name);
+
+ %% postprocess
+ % store result in structure for saving now
+ % only apply postprocessing to final pass
+ fprintf(' postprocessing...\n');
fs = struct();
- % only save final pass
for c = 1 : n_correl
vname = sprintf('C%d',c);
- VC = piv_result{c}(end);
+ VC = scaledt(piv_result{c}(end), piv_dt);
+ VC = vector_postprocess(VC, operations);
VC.attributes= atts;
fs.(vname) = VC;
end
+
+ %% save result
+ t_save = tic;
+ fprintf(' save to %s\n', vec_name);
if exist(vec_name, 'file')
warning('%s already exists, overwriting result', vec_name);
end
@@ -191,19 +221,24 @@ function job_execute(job_file)
%% report
t_samp = toc(t_samp);
- fprintf(' took %0.1f s [%0.1f load, %0.1f preproc, %0.1f proc, %0.1f save]\n', t_samp, t_load, t_preproc, t_proc, t_save);
+ fprintf(' took %0.2f s [%0.2f load, %0.2f preproc, %0.2f proc, %0.2f save]\n', t_samp, t_load, t_preproc, t_proc, t_save);
fprintf(' %2d h %2d m %2d s elapsed so far\n', hms(toc(t_begin)));
+ %% figure settings
+ offset = 10;
+
+
%% plot vector field
pass = length(piv_result{1});
if plt_opts.b_plot_u
- figure(1);
+ figure(1+offset);
clf;
ax = [];
for c = 1 : n_correl
V = piv_result{c}(pass);
+ V.ux(V.nan_mask)= nan;
subplot(1,n_correl,c);
- ax(c) = gca;
+ ax(end+1) = gca;
imagesc(V.win_x, V.win_y, V.ux')
axis equal tight ij;
colormap(redblue);
@@ -217,18 +252,19 @@ function job_execute(job_file)
%% plot v component field
if plt_opts.b_plot_v
- figure(2);
+ figure(2+offset);
clf;
ax = [];
for c = 1 : n_correl
V = piv_result{c}(pass);
+ V.uy(V.nan_mask)= nan;
subplot(1,n_correl,c);
- ax(c) = gca;
+ ax(end+1) = gca;
imagesc(V.win_x, V.win_y, V.uy')
axis equal tight ij;
colormap(redblue);
colorbar;
- set(gca,'clim',[-4 4]);
+ set(gca,'clim',[-16 16]);
title(sprintf('v component, camera %d', c));
end
linkaxes(ax,'xy');
@@ -237,13 +273,14 @@ function job_execute(job_file)
%% plot correlation value
if plt_opts.b_plot_cv
- figure(4);
+ figure(4+offset);
clf;
ax = [];
for c = 1 : n_correl
V = piv_result{c}(pass);
subplot(1,n_correl,c);
- ax(c) = gca;
+ ax(end+1) = gca;
+ V.peak_mag(V.nan_mask)= nan;
imagesc(V.win_x, V.win_y, V.peak_mag')
axis equal tight ij;
hold on;
@@ -258,6 +295,7 @@ function job_execute(job_file)
end
linkaxes(ax,'xy');
drawnow;
+ fprintf('mean correlation value = %0.2f\n', nanmean(V.peak_mag(:)));
end
end
end
\ No newline at end of file
diff --git a/job_postprocess.m b/job_postprocess.m
new file mode 100644
index 0000000000000000000000000000000000000000..3fdd3555f38035a946aae5e46b4c15d023e40fbf
--- /dev/null
+++ b/job_postprocess.m
@@ -0,0 +1,65 @@
+function job_postprocess(job_file)
+ % load a PIV processing job and execute it
+
+ %% load options from job file
+ job = importdata(job_file);
+ % source images
+ src_dir = job.src_dir;
+ src_pat = job.src_pat;
+ src_range = job.src_range;
+ n_src = length(src_range);
+ % output vector fields
+ b_overwrite = job.b_overwrite;
+ dst_dir = job.dst_dir;
+ dst_pat = job.dst_pat;
+ % postprocessing operations to perform
+ operations = job.operations;
+
+ %% report
+ fprintf('--------------------------------------------------------------\n');
+ fprintf('PROCESSING START %s\n', datestr(now));
+ t_begin = tic;
+
+
+ %% loop over files to be processed
+ for n = 1 : n_src
+ %% load
+ t_samp = tic;
+ % check source exists
+ src_name = sprintf(src_pat, src_dir, src_range(n));
+ fprintf('[%3d/%3d] %s \n', n, n_src, src_name);
+ if ~exist(src_name, 'file')
+ warning('source %s does not exist, skipping\n', src_name);
+ continue;
+ end
+ % check output does not exist
+ dst_name = sprintf(dst_pat, dst_dir, src_range(n));
+ if exist(dst_name, 'file') && ~b_overwrite
+ warning('output %s already exists, skipping\n', dst_name);
+ continue;
+ end
+ % load vector field from file
+ t_load = tic;
+ V = load(src_name);
+ t_load = toc(t_load);
+
+ %% apply postprocessing chain to vector field
+ t_proc = tic;
+ W = vector_postprocess(V.C1, operations);
+ t_proc = toc(t_proc);
+
+ %% save result
+ t_save = tic;
+ fprintf(' save to %s\n', dst_name);
+ if exist(dst_name, 'file')
+ warning('%s already exists, overwriting result', dst_name);
+ end
+ save(dst_name, '-struct', 'W');
+ t_save = toc(t_save);
+
+ %% report
+ t_samp = toc(t_samp);
+ fprintf(' took %0.1f s [%0.1f load, %0.1f proc, %0.1f save]\n', t_samp, t_load, t_proc, t_save);
+ fprintf(' %2d h %2d m %2d s elapsed so far\n', hms(toc(t_begin)));
+ end
+end
\ No newline at end of file
diff --git a/misc/gradient_o2_2d.m b/misc/gradient_o2_2d.m
new file mode 100644
index 0000000000000000000000000000000000000000..b1f07b27ca369291ef4f96df9c17f111aa482f51
--- /dev/null
+++ b/misc/gradient_o2_2d.m
@@ -0,0 +1,33 @@
+function [fx, fy] = gradient_o4_2d(f, dx, dy)
+ % calculate second order accurate approximation to gradient of f(x,y)
+ % using the formula
+ % df/dx = (u(i+1,j,k) - u(i-1,j,k))/(2dx)
+
+ %% fill in missing arguments
+ if nargin < 2
+ dx = 1;
+ end
+ if nargin < 3
+ dy = dx;
+ end
+
+ %% generate kernels
+ kernel_x = reshape([1 0 -1] / (2*dx), [3 1]);
+ kernel_y = reshape([1 0 -1] / (2*dy), [1 3]);
+
+ %% convolve over given dimension
+ fx = convn(f, kernel_x, 'same');
+ fy = convn(f, kernel_y, 'same');
+
+ %% mask invalid regions with nan
+ if size(f,1) >= 3
+ fx([1 end],:) = nan;
+ else
+ fx = nan(size(f));
+ end
+ if size(f,2) >= 3
+ fy(:,[1 end]) = nan;
+ else
+ fy = nan(size(f));
+ end
+end
\ No newline at end of file
diff --git a/misc/gradient_o4_2d.m b/misc/gradient_o4_2d.m
new file mode 100644
index 0000000000000000000000000000000000000000..9f3844765c6ed9908aea848bb855dfa6110621c8
--- /dev/null
+++ b/misc/gradient_o4_2d.m
@@ -0,0 +1,34 @@
+function [fx, fy] = gradient_o4_2d(f, dx, dy)
+ % calculate fourth order accurate approximation to gradient of f(x,y)
+ % using the formula
+ % df/dx = (-u(i+2,j,k) + 8u(i+1,j,k) - 8u(i-1,j,k) + u(i-2,j,k))/(12dx)
+ % df/dy = (-u(i,j+2,k) + 8u(i,j+1,k) - 8u(i,j-1,k) + u(i,j-2,k))/(12dy)
+
+ %% fill in missing arguments
+ if nargin < 2
+ dx = 1;
+ end
+ if nargin < 3
+ dy = dx;
+ end
+
+ %% generate kernels
+ kernel_x = reshape([-1 8 0 -8 1] / (12*dx), [5 1 1]);
+ kernel_y = reshape([-1 8 0 -8 1] / (12*dy), [1 5 1]);
+
+ %% convolve over given dimension
+ fx = convn(f, kernel_x, 'same');
+ fy = convn(f, kernel_y, 'same');
+
+ %% mask invalid regions with nan
+ if size(f,1) > 5
+ fx([1:2 end-1 end],:)= nan;
+ else
+ fx = nan(size(f));
+ end
+ if size(f,2) > 5
+ fy(:,[1:2 end-1 end])= nan;
+ else
+ fy = nan(size(f));
+ end
+end
\ No newline at end of file
diff --git a/misc/readimg.m b/misc/readimg.m
new file mode 100644
index 0000000000000000000000000000000000000000..093ab159237d4639f4eaac11fd974f5786ab5b4a
--- /dev/null
+++ b/misc/readimg.m
@@ -0,0 +1,56 @@
+function [im_ary, atts] = readimg(file_name)
+ [~,~,ext] = fileparts(file_name);
+ if strcmpi(ext, '.imx')
+ %% load imx image array
+ imx = readimx(src_name);
+ im_ary = cellfun(@(f) single(f.Components{1}.Planes{1}), imx.Frames, 'UniformOutput', false);
+ atts = readAttributes(imx);
+ elseif strcmpi(ext, '.tif') || strcmpi(ext, '.tiff')
+ %% load TIF image array
+ info = imfinfo(file_name);
+ n_pages = length(info);
+ im_ary = cell(1, n_pages);
+ for n = 1 : n_pages
+ im_ary{n} = single(imread(file_name, n));
+ end
+
+ %% read optional attributes from TIFF XMP tag
+ % using custom headers from gige data logger
+ t = Tiff(file_name, 'r');
+ xmlstr = '<metadata></metadata>';
+ try
+ xmlstr = sprintf('<metadata>\n%s\n</metadata>', t.getTag('XMP'));
+ end
+ dom = xmlstr2struct(xmlstr);
+ metadata = dom.metadata;
+ atts = containers.Map;
+ % parse fields, if present
+ if isfield(metadata, 'AcqTime')
+ atts('AcqTime') = metadata.AcqTime.Text;
+ end
+ if isfield(metadata, 'DeltaT')
+ % hard coded to use microseconds as units for now
+ % unit information is provided in
+ units_to_seconds= 1E-6;
+ atts('DeltaT') = str2double(metadata.DeltaT.Text) * units_to_seconds;
+ end
+ if isfield(metadata, 'Exposure')
+ units_to_seconds= 1E-6;
+ atts('Exposure')= str2double(metadata.Exposure.Text) * units_to_seconds;
+ end
+ if isfield(metadata, 'SampleIndex')
+ atts('SampleIndex') = str2double(metadata.SampleIndex.Text);
+ end
+ elseif strcmpi(ext, '.jpg')
+ %% load single JPG image
+ img = imread(file_name);
+ % convert to grayscale if necessary
+ img = sqrt(sum(img.^2, 3));
+ im_ary = {img};
+ atts = containers.Map;
+ atts('SampleIndex')= 1;
+ atts('DeltaT') = 1;
+ else
+ error('unrecognised file format');
+ end
+end
\ No newline at end of file
diff --git a/misc/scale2d.m b/misc/scale2d.m
new file mode 100644
index 0000000000000000000000000000000000000000..70b40b79d782d22669bd788a086ff52c1abf02b8
--- /dev/null
+++ b/misc/scale2d.m
@@ -0,0 +1,64 @@
+function new = scale2d(old, calib, b_keep)
+ % W = scale2d(old, calib)
+ % rescale a 2D vector field of pixel displacement by a specified calibration
+ %
+ % Input:
+ % old vector field structure with image space coordinates
+ % calib calibration structure
+ % b_keep flag to keep pixel displacement information
+ %
+ % Output:
+ % new vector field with components rescaled to object space
+ %
+ % calibration structure contains:
+ % model char (optional) string describing type of calibration model
+ % P 3x3 linear transformation matrix
+ %
+ % currently, only a linear transformation is implemented
+ % see p2d_im2obj and p2d_obj2im
+ %
+ % linear transform is defined
+ % lambda * [x;y;1] = P*[X;Y;1]
+ % where
+ % (x,y) coordinates in image space
+ % (X,Y) coordinates in object space
+ % lambda scale factor
+
+ %% extract velocity in image coordinates
+ im_u = old.ux;
+ im_v = old.uy;
+
+ %% apply transformation to position
+ [im_x, im_y]= ndgrid(old.win_x, old.win_y);
+ n_vec = numel(im_x);
+ n_windows = size(im_x);
+ im_pos = [im_x(:) im_y(:)]';
+ gl_pos = p2d_map_im2obj(calib.P, im_pos);
+ X = reshape(gl_pos(1,:), n_windows);
+ Y = reshape(gl_pos(2,:), n_windows);
+
+ %% apply transformation to displacement
+ im_aug = [im_pos; ones(1, n_vec)];
+ UV_aug = zeros(3, n_vec);
+ P = calib.P;
+ lambda = P(3,:)*[gl_pos; ones(1, n_vec)];
+ for k = 1 : n_vec
+ P_aug = [P(:,1:2), -im_aug(:,k)];
+ UV_aug(:,k) = lambda(k) * (P_aug \ [im_u(k); im_v(k); 0]);
+ end
+ U = reshape(UV_aug(1,:), n_windows);
+ V = reshape(UV_aug(2,:), n_windows);
+
+ %% update fields
+ % discard old fields
+ new = old;
+ if nargin < 3 || ~b_keep
+ new = rmfield(new, {'ux', 'uy', 'win_x', 'win_y'});
+ end
+ % add new fields
+ new.X = X;
+ new.Y = Y;
+ new.U = U;
+ new.V = V;
+ new.calib = calib;
+end
\ No newline at end of file
diff --git a/misc/scaledt.m b/misc/scaledt.m
new file mode 100644
index 0000000000000000000000000000000000000000..3a9c484e727e3067588877577dadc4b242a465f5
--- /dev/null
+++ b/misc/scaledt.m
@@ -0,0 +1,28 @@
+function new = scaledt(old, new_dt)
+ % W = scaledt(V, scale)
+ % rescale PIV delta-t by new value
+ % U -> U * dt / new_dt
+ % dt -> new_dt
+
+ %% define new scale factor
+ new = old;
+ scale = old.dt / new_dt;
+ new.dt = new_dt;
+
+ %% apply scale to velocity components if present
+ if isfield(old, 'ux')
+ new.ux = old.ux * scale;
+ end
+ if isfield(old, 'uy')
+ new.uy = old.uy * scale;
+ end
+ if isfield(old, 'U')
+ new.U = old.U * scale;
+ end
+ if isfield(old, 'V')
+ new.V = old.V * scale;
+ end
+ if isfield(old, 'W')
+ new.W = old.W * scale;
+ end
+end
\ No newline at end of file
diff --git a/misc/test_scale2d.m b/misc/test_scale2d.m
new file mode 100644
index 0000000000000000000000000000000000000000..85ac018c362bac52d0a16a681483a7fc48b3e968
--- /dev/null
+++ b/misc/test_scale2d.m
@@ -0,0 +1,42 @@
+
+clear;
+load('F:\mixer\231020\calibration\CAM1-01.mat');
+
+%% create a random velocity field
+win_x = 16 : 16 : 6600;
+win_y = 16 : 16 : 4400;
+[im_x,im_y] = ndgrid(win_x, win_y);
+im_pos = [im_x(:) im_y(:)]';
+obj_pos = p2d_map_im2obj(P, im_pos);
+n_windows = size(im_x);
+
+% rando velocity
+im_ux = randn(n_windows);
+im_uy = randn(n_windows);
+% construct object space velocity using numerical approximation
+im_pos_plus = im_pos + [im_ux(:) im_uy(:)]'/2;
+im_pos_minus= im_pos - [im_ux(:) im_uy(:)]'/2;
+obj_vel = p2d_map_im2obj(P, im_pos_plus) - p2d_map_im2obj(P, im_pos_minus);
+obj_ux = reshape(obj_vel(1,:), n_windows);
+obj_uy = reshape(obj_vel(2,:), n_windows);
+
+%% apply rescaling
+V_im = struct('win_x', win_x, 'win_y', win_y, ...
+ 'ux', im_ux, 'uy', im_uy);
+calib = struct('P', P);
+V_obj = scale2d(V_im, calib, false);
+
+%% compare
+rms(obj_ux(:) - V_obj.U(:))
+rms(obj_uy(:) - V_obj.V(:))
+
+%% compare
+figure(1);
+clf;
+plot(obj_ux(:), V_obj.U(:), '.');
+hold on;
+plot(obj_uy(:), V_obj.V(:), '.');
+axis equal tight;
+
+
+
diff --git a/misc/xmlstr2struct.m b/misc/xmlstr2struct.m
new file mode 100644
index 0000000000000000000000000000000000000000..79aed7eb765499b15c9ed168fc7ac8f29a04a957
--- /dev/null
+++ b/misc/xmlstr2struct.m
@@ -0,0 +1,28 @@
+function s = xmlstr2struct(xml)
+ % s = xmlstr2struct(xml)
+ % parse an XML string into a MATLAB structure
+ %
+ % Input:
+ % xml xml markup string
+ % Output:
+ % s MATLAB structure containing parsed output of xml markup
+
+ % create a temporary file name
+ s = [];
+ fname = tempname;
+ try
+ %% create temporary file
+ fp = fopen(fname, 'w');
+ fwrite(fp, xml);
+ fclose(fp);
+
+ %% read xml
+ s = xml2struct(fname);
+
+ %% delete
+ delete(fname);
+ catch
+ % delete temporary file
+ delete(fname);
+ end
+end
\ No newline at end of file
diff --git a/sampledata/NACA0020_00000.jpg b/sampledata/NACA0020_00000.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..aeab5427f7ee8eedb512cb77d94b4a2c29ac7be0
Binary files /dev/null and b/sampledata/NACA0020_00000.jpg differ
diff --git a/sampledata/NACA0020_00001.jpg b/sampledata/NACA0020_00001.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..43be1936f31e5d87212f660a2c5be0275ac61d00
Binary files /dev/null and b/sampledata/NACA0020_00001.jpg differ
diff --git a/vector_postprocess.m b/vector_postprocess.m
new file mode 100644
index 0000000000000000000000000000000000000000..b26d7dcf10c5939a80f3257a9fd28374af933332
--- /dev/null
+++ b/vector_postprocess.m
@@ -0,0 +1,44 @@
+function W = vector_postprocess(V, ops)
+ % apply postprocessing to vector field
+ % by sequentially applying operations in ops to V
+ %
+ % each element in ops has at least one member, 'type', which determines what
+ % type of filter is to be applied. it is a string which selects from a
+ % number of different available filters
+ %
+ % Type/member Default Description
+ % ----------------------------------------------------------------------
+ % scale2d apply scaling from image to object space
+ % .P eye(3) 2D projection matrix representing mapping
+ % from image space to object space
+ % .keep false flag to indicate whether old pixel
+ % displacement data should be retained
+ % scaledt rescale PIV delta t
+ % .dt 1 new value of PIV delta-t to rescale data to
+
+ %% iteratively apply filters
+ W = V;
+ n_ops = length(ops);
+
+ for i = 1 : n_ops
+ operation = ops{i};
+ if strcmp(operation.type, 'scale2d')
+ %% scale from image to object space
+ calib = struct('P', eye(3));
+ keep = false;
+ if isfield(operation, 'P'), calib.P = operation.P; end
+ if isfield(operation, 'keep'), keep = operation.keep; end
+ W = scale2d(W, calib, keep);
+ elseif strcmp(operation.type, 'scaledt')
+ %% scale to new delta t
+ dt = 1; % scales to pixel displacement (unit delta t)
+ if isfield(operation, 'dt'), dt = operation.dt; end
+ W = scaledt(W, dt);
+ elseif strcmp(operation.type, 'null')
+ %% null operation, do nothing
+ else
+ error('unrecognised filter type %s', filter.type);
+ end
+ end
+end
+
\ No newline at end of file