diff --git a/image_preprocess.m b/image_preprocess.m index d9d762cfec9a31536ccd043483a4b316a7fbc289..42fe1a2650580a88f161f74469abb0c4c4ea0abe 100644 --- a/image_preprocess.m +++ b/image_preprocess.m @@ -70,6 +70,9 @@ function im_filt = image_preprocess(im, filters) if isfield(filter, 'n'), n = filter.n; end thresh = im_med + n*im_std; im_filt = min(im_filt, thresh); + elseif strcmp(filter.type, 'null') + %% null filter, do nothing + im_filt = im_filt; else error('unrecognised filter type %s', filter.type); end diff --git a/job_create.m b/job_create.m index 579a157fbe9244a38417ad50466ecfaa153d5060..199fccb6301079d25c2e263593be21954640e0a6 100644 --- a/job_create.m +++ b/job_create.m @@ -9,31 +9,35 @@ addpath(genpath(fileparts(mfilename('fullpath')))); %% input parameters: file paths % where to save jobs -%job_idx = 1; -job_dir = 'H:\davis\multiscale\Gen1_12\jobs\'; +job_idx = 1; +job_dir = 'J:\FEE\EFMRL\John\TR-IM7\jobs\'; 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 -%src_dir = 'H:\davis\multiscale\Gen1_12\U10\U10_x-axis=257.988_x-axis=257.987_x-axis=229.986\'; +src_dir = 'J:\FEE\EFMRL\John\TR-IM7'; src_pat = '%s/B%05d.im7'; % source image to process -src_range = 1:1000; % range of source images to process +src_range = 1:3; % 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) msk_dir = [src_dir '/MASK/']; % mask image directory -msk_pat = '%s/B0001.im7'; % pattern for mask image input (ignore if not using) +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 '/C12']; % directory to save vector fields to vec_pat = '%s/B%05d.mat'; % output vector field file pattern %% input parameters: processing % processing flags -b_use_mask = true; % use mask file? +b_use_mask = false; % use mask file? b_overwrite = true; % overwrite existing vector fields? -b_execute_now = false; % start processing now? +b_execute_now = true; % start processing now? % PIV correlation settings -i_frames = {[1 2],[3 4]}; % indices of frame pairs to cross-correlate +i_frames = {[1 1],[2 2], [3 3]}; % 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 @@ -52,27 +56,27 @@ min_filt_size = [17 17]; % sliding bg filter size gauss_filt_size = [3 3]; % gaussian filter size norm_filt_size = [17 17]; % normalisation filter size -filters = {}; +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 Gaussian filter -filters{end+1} = struct('type', 'gaussian', 'size', gauss_filt_size, 'sigma', 0.67); -% min-max normalisation filter -filters{end+1} = struct('type', 'norm', 'size', norm_filt_size, 'max_gain', 1/32); +% % 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 Gaussian filter +% filters{end+1} = struct('type', 'gaussian', 'size', gauss_filt_size, 'sigma', 0.67); +% % min-max normalisation filter +% filters{end+1} = struct('type', 'norm', 'size', norm_filt_size, 'max_gain', 1/32); %% 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 = false; % plot source images? -plt_opts.b_plot_u = false; % plot u component field? -plt_opts.b_plot_v = false; % plot v component field? -plt_opts.b_plot_cv = false; % plot correlation value? +plt_opts.font_size = 24; % figure font size +plt_opts.b_plot_src = true; % plot source images? +plt_opts.b_plot_u = true; % plot u component field? +plt_opts.b_plot_v = true; % plot v component field? +plt_opts.b_plot_cv = true; % plot correlation value? %% check job directory job_completed_dir = [job_dir '/completed/']; @@ -101,6 +105,11 @@ if exist(job_file, 'file') 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); @@ -187,6 +196,8 @@ job = struct(); 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; diff --git a/job_execute.m b/job_execute.m index 3c25f31b8f117011fa1be148e2e12fb694f364a9..6d1d3d6f2718fc7f27f2b19d6be11e812561b8af 100644 --- a/job_execute.m +++ b/job_execute.m @@ -8,6 +8,8 @@ function job_execute(job_file) src_pat = job.src_pat; src_range = job.src_range; n_src = length(src_range); + % timeseries PIV flag + b_timeseries = job.b_timeseries; % background image b_use_bg = job.b_use_bg; bg_dir = job.bg_dir; @@ -89,8 +91,13 @@ function job_execute(job_file) % read attributes atts = readAttributes(imx); - s = atts('Reference time dt 1'); - piv_dt = 1E-6 * s.value; + 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 t_load = toc(t_samp); %% preprocess images @@ -109,17 +116,45 @@ function job_execute(job_file) end t_preproc = toc(t_preproc); - - %% PIV cross-correlation + %% create image pairs + im_pairs = cell(2, n_correl); + im_correl_mask = cell(1, n_correl); + if b_timeseries + %% time-series cross-correlation + % preload first image + if n == 1 + im_filt_prev= im_filt; + continue; + end + % allocate images from timeseries to image pairs + for c = 1 : n_correl + iA = piv_opts(c).i_frames(1); + iB = piv_opts(c).i_frames(2); + im_pairs(1,c) = im_filt_prev(iA); + im_pairs(2,c) = im_filt(iB); + im_correl_mask{c} = im_mask_ary{iA} | im_mask_ary{iB}; + end + % save current pair for next iteration + im_filt_prev = im_filt; + else + %% regular PIV + for c = 1 : n_correl + iA = piv_opts(c).i_frames(1); + iB = piv_opts(c).i_frames(2); + im_pairs(1,c) = im_filt(iA); + im_pairs(2,c) = im_filt(iB); + im_correl_mask{c} = im_mask_ary{iA} | im_mask_ary{iB}; + end + end + + %% cross-corrleate piv_result = cell(1,n_correl); t_proc = tic; for c = 1 : n_correl %% get images to cross-correlate - iA = piv_opts(c).i_frames(1); - iB = piv_opts(c).i_frames(2); - im_A = im_filt{iA}; - im_B = im_filt{iB}; - + im_A = im_pairs{1,c}; + im_B = im_pairs{2,c}; + %% plot source, if requested if plt_opts.b_plot_src figure(10); @@ -128,11 +163,11 @@ function job_execute(job_file) axis equal tight ij; drawnow; end - + %% cross-correlate fprintf(' correlation %d of %d\n', c, n_correl); piv_result{c} = PIV_2D_wdef(im_A, im_B, ... - im_mask_ary{iA} | im_mask_ary{iB}, ... + im_correl_mask{c}, ... piv_opts(c)); end t_proc = toc(t_proc); @@ -211,6 +246,11 @@ function job_execute(job_file) ax(c) = gca; imagesc(V.win_x, V.win_y, V.peak_mag') axis equal tight ij; + hold on; + [xm,ym] = ndgrid(V.win_x, V.win_y); + ir = 1 : 4 : V.n_windows(1); + jr = 1 : 4 : V.n_windows(2); + quiversc(xm(ir,jr),ym(ir,jr),V.ux(ir,jr),V.uy(ir,jr),10,'k'); colormap(jet); colorbar; set(gca,'clim',[0 1]);