From 7959c1c563453d11b3c4f31ab38eb63a9fa3d026 Mon Sep 17 00:00:00 2001 From: "Jun, Brian H" Date: Tue, 27 Apr 2021 21:24:00 -0400 Subject: [PATCH] Add files via upload --- Hessian_2D.m | 53 ++++ ParabolaRoots.m | 57 ++++ calculate_GSM.m | 364 +++++++++++++++++++++++++ cell_connect_tracks.m | 134 +++++++++ cell_core_segmentation.m | 223 +++++++++++++++ cell_downsample_Hessian_segmentation.m | 158 +++++++++++ cell_parameterization.m | 166 +++++++++++ cell_plot_sequence.m | 78 ++++++ cell_postprocessing.m | 271 ++++++++++++++++++ cell_processing_main.m | 136 +++++++++ cell_refinement.m | 129 +++++++++ cell_tracking.m | 32 +++ distinguishable_colors.m | 152 +++++++++++ eig2.m | 34 +++ ellipse.m | 149 ++++++++++ finding_true_num_2D.m | 19 ++ parsave.m | 3 + polartrans.m | 89 ++++++ socdiff.m | 56 ++++ weighted_nearest_neighbor3D_MP.m | 137 ++++++++++ 20 files changed, 2440 insertions(+) create mode 100644 Hessian_2D.m create mode 100644 ParabolaRoots.m create mode 100644 calculate_GSM.m create mode 100644 cell_connect_tracks.m create mode 100644 cell_core_segmentation.m create mode 100644 cell_downsample_Hessian_segmentation.m create mode 100644 cell_parameterization.m create mode 100644 cell_plot_sequence.m create mode 100644 cell_postprocessing.m create mode 100644 cell_processing_main.m create mode 100644 cell_refinement.m create mode 100644 cell_tracking.m create mode 100644 distinguishable_colors.m create mode 100644 eig2.m create mode 100644 ellipse.m create mode 100644 finding_true_num_2D.m create mode 100644 parsave.m create mode 100644 polartrans.m create mode 100644 socdiff.m create mode 100644 weighted_nearest_neighbor3D_MP.m diff --git a/Hessian_2D.m b/Hessian_2D.m new file mode 100644 index 0000000..3d05710 --- /dev/null +++ b/Hessian_2D.m @@ -0,0 +1,53 @@ +function [lambda1,lambda2] = Hessian_2D(image) + +% ========================================================================= +% sub function called by iterative particle reconstruction +% calculate the Hessian maxtrix of a 2D image +% and its eigenvalues for filtering. +% by Tianqi Guo, Aug 2017 +% ========================================================================= + +% display('Hessian filtering') + +% addpath('Eig3Folder/'); + +% Gaussian smoothing parameters for eigenvalue matrix +filter_sd = 0.65; + +% first order derivatives of intensity volume +gx = socdiff(image,[],2); +gy = socdiff(image,[],1); + +% second order derivatives of intensity volume +gxx = socdiff(gx,[],2); +gxy = socdiff(gx,[],1); + +gyx = socdiff(gy,[],2); +gyy = socdiff(gy,[],1); + +% eigenvalues of the Hessian matrix at each voxel +% sorted as absolute values lambda1 < lambda2 + +H = cat(4,gxx,gxy,gyx,gyy); + +H = reshape(H,[size(image(:),1),2,2]); + +H = permute(H,[3 2 1]); + +% eigenvalue for Hessian matrix and its absolute value +e = eig2(H); +e_abs = abs(e); + +[~, ind] = sort(e_abs); +[m,n] = size(e); +e_sort = real(e(bsxfun(@plus,ind,(0:n-1)*m))); + +lambda1 = reshape(e_sort(1,:),size(image)); +lambda2 = reshape(e_sort(2,:),size(image)); + +% Gaussian smoothing of eigenvalue matrices to suppress noise +lambda1 = imgaussfilt(lambda1,filter_sd); +lambda2 = imgaussfilt(lambda2,filter_sd); + + +end \ No newline at end of file diff --git a/ParabolaRoots.m b/ParabolaRoots.m new file mode 100644 index 0000000..6a2d559 --- /dev/null +++ b/ParabolaRoots.m @@ -0,0 +1,57 @@ +function roots = ParabolaRoots(varargin) +% function roots = ParabolaRoots(P) +% +% Find roots of second order polynomials +% +% INPUT +% P: (n x 2) array, each row corresponds to coefficients of each +% polynomial, P(:,1)*x^2 + P(:,2)*x + P(:,3) +% OUTPUT +% roots: (n x 2) array, each row correspond to the roots of P +% +% To adjust the parameter below which the the discriminant is considerered +% as nil, use +% ParabolaRoots(P, tol) +% Adjusting tol is useful to avoid the real roots become complex due to +% numerical accuracy. The default TOL is 0 +% +% See also: roots, CardanRoots, eig2 +% +% Author: Bruno Luong +% History: +% Original 27-May-2010 + +% Adjustable parameter +tol = 0; + +if nargin<3 + P = varargin{1}; + a = P(:,1); + b = P(:,2); + c = P(:,3); + if nargin>=2 + tol = varargin{2}; + end +else + [a b c] = deal(varargin{1:3}); + if nargin>=4 + tol = varargin{2}; + end +end + +if ~isequal(a,1) + b = b./a; + c = c./a; +end + +b = 0.5*b; + +delta = b.^2 - c; +delta(abs(delta) 1 + t = [s(:,1) s(:,2) s(:,1)+s(:,3) s(:,2)+s(:,4)]; %[left bottom right top] + leftBottomCorner = min(t(:,1:2)); + rightTopCorner = max(t(:,3:4)); + MBB_xy = [leftBottomCorner rightTopCorner-leftBottomCorner]; %[left bottom width height] +else + MBB_xy = s; +end + +rectangularity = sum(area_rotated) / (MBB_xy(3)*MBB_xy(4)); + +%========================================================================== +% Polar transform based + +% Polar transform +nrad = 300; +ntheta = 360; +linlog = 'linear'; +shape = 'valid'; +[theta, radius, pim] = polartrans(BW_img, nrad, ntheta, CoM(1), CoM(2), linlog, shape); + +% r(thetha) + +r_theta = zeros(1,ntheta); +for t = 1:ntheta + %ind = find(pim(:,t)==0,1,'first'); + ind = find(pim(:,t)==1,1,'last'); + if isempty(ind) + r_theta(t) = 0; + else + r_theta(t) = radius(ind(1),t); + end +end +r_theta = smooth(r_theta(:)); +r_avg = mean(r_theta); +r_std = var(r_theta).^0.5; +r_max = max(r_theta); + +% circularities +circularity_Ravg = 1 - trapz(theta(1,:),(r_theta - r_avg).^2)/(r_avg^2*trapz(theta(1,:),ones(size(theta(1,:))))); +circularity_Rmax = area/(pi*r_max^2); +circularity_Rstd = 1 - r_std/r_avg; + +% ========================================================================= +% if orient_degree < 0 +% orient_degree = orient_degree + 90; +% end +measures = [IWC(1), IWC(2), CoM(1), CoM(2), I_avg, I_max... + area, per, ra, rb, r_eq, orient_degree,... + elongation, eccentricity_chord, eccentricity_ellipse, elongatedness, aspectRatio, ... + circularity_AP, circularity_Rstd, circularity_Ravg, circularity_Rmax, ... + ellipticity, squareness, triangularity, rectangularity, convexity_a, convexity_p]; +cats = {'Basics','IWC','Center of mass','I_mean','I_max','Area (pixel^2)', 'Perimeter (pixel)',... + 'Major axis (pixel)', 'Minor axis (pixel)','Equiv R (pixel)', 'Orientation (°)', ... + 'Eccentricities','Elongation (FFT)','Ecc (chord)','Ecc (ellipse)', 'Elongatedeness', 'Ratio (Rb/Ra)',... + 'Circularities','Circ (AP)', 'Circ (r_std)', 'Circ (r_avg)', 'Circ (r_max)', ... + 'Shape factors','Ellipticity','Squareness', 'Triangularity', 'Rectangularity', 'Convexity (A)', 'Convexity (P)'}; +table_data = {}; +table_data{1} = []; +table_data{2} = sprintf('(%.1f,%.1f)',IWC); +table_data{3} = sprintf('(%.1f,%.1f)',CoM); +for i = 4 : 11 + table_data{i} = sprintf('%.1f',measures(i+1)); +end +table_data{12} = []; +for i = 13 : 17 + table_data{i} = sprintf('%.4f',measures(i)); +end +table_data{18} = []; +for i = 19 : 22 + table_data{i} = sprintf('%.4f',measures(i-1)); +end +table_data{23} = []; +for i = 24 : 29 + table_data{i} = sprintf('%.4f',measures(i-2)); +end +if optionDebug + file_list = dir(fullfile('debug','*.tif')); + N = length(file_list) + 1; + + close all; + fig = figure('position',[00 200 1500 900]); + + % subplot(241) + % imshow(BW_rotated); hold on; + % rectangle('Position',MBB_xy,'edgecolor','m','linestyle',':','linewidth',2); + % title('Minimum bounding rectangle') + % + % subplot(242) + % imshow(BW_img_rotated); hold on; + % plot(chords_xy(:,1),chords_xy(:,2),'-m','linewidth',2); + % plot(chords_xy(:,3),chords_xy(:,4),'-m','linewidth',2) + % title('Longest chord') + + subplot(243) + s = size(raw_img); + imagesc(uint8(raw_img)); axis equal off; hold on; + h1 = plot(CoM(1),CoM(2),'bO','markersize',5,'linewidth',2); + h2 = ellipse(rb,ra,pi/2-orient_radian,CoM(1),CoM(2),'r',50); hold on; + h3 = plot(Con_hul(:,1),Con_hul(:,2),'g-.','linewidth',2); + h4 = plot([CoM(1) s(2)],[CoM(2) CoM(2)],'m--','linewidth',2); hold on + plot([CoM(1) s(2)],[CoM(2) s(1)- (CoM(2)+(s(2)-CoM(1))*tan(orient_radian))],'m--','linewidth',2); + title('Segmented image') + legend([h1,h2,h3,h4],... + {'Center of mass','Equivalent eclipse','Smallest enclosing convex','Orientation angle'},... + 'location','southeast') + xlim([1 s(2)]) + ylim([1 s(1)]) + + subplot(245) + contourf(theta,radius,pim,'edgecolor','none'); hold on; + plot(theta(1,:),r_theta,'r-','linewidth',2) + xlabel('\theta (Radian)') + ylabel('r (pixel)') + title('Polar transform') + + subplot(246) + plot(boundary_x,'.'); hold on; + plot(boundary_y,'.'); hold on; + grid on + title('Boundary coordinates') + legend('x','y','location','south') + xlabel('Number of points') + ylabel('Coordinates (pixel)') + xlim([1,n_boundary]) + + subplot(247) + n_2 = floor(length(c)/2); + if mod(n_boundary,2) == 0 + n = -n_2:(n_2-1); + else + n = -n_2:n_2; + end + r_0 = (x_0^2 + y_0^2)^0.5; + stem(n,fftshift(c),'filled','linewidth',2,'markersize',5) + xlim([-10 10]) + ylim([0 r_0*1.1]) + grid on; + xlabel('Descriptor number (n)') + ylabel('|a_n + i * b_n|') + title('Complex Fourier transform') + + subplot(2,4,[4,8]) + %cnames = {'Measures','Values'}; + t = uitable(fig,'Data',cat(2,cats(:),table_data(:)),... + 'ColumnName',[],'RowName',[],'ColumnWidth',{120,100},... + 'ColumnFormat',{'char','char'}); + t.FontSize = 10; + plot(3) + pos = get(subplot(2,4,[4,8]),'position'); + delete(subplot(2,4,[4,8])) + set(t,'units','normalized') + set(t,'position',pos) + % jScroll = findjobj(t); + % jTable = jScroll.getViewport.getView; + % jTable.setAutoResizeMode(jTable.AUTO_RESIZE_SUBSEQUENT_COLUMNS); + % drawnow; + + % print(fullfile('debug',sprintf('%0.4d.tif',N)),'-dtiff','-r100') + +end +end + +function M = raw_moments(BW_img,p,q) +s = size(BW_img); +[X,Y] = meshgrid(1:s(2),1:s(1)); +integrand = X.^p .* Y.^q .* BW_img; +M = sum(integrand(:)); +end + +function [BW_img_rotated,chords_xy] = longest_chord(BW_img,boundary_x,boundary_y) +%========================================================================== +% try +s = size(BW_img); +x0 = round(s(2)/2); +y0 = round(s(1)/2); +n_boundary = length(boundary_x); + +chord_length_max = -1; +for i = 1:n_boundary-1 + for j = i+1:n_boundary + flag = 0; + if boundary_x(i) == boundary_x(j) + step = sign(round(boundary_y(j))-round(boundary_y(i))); + for yy = round(boundary_y(i)) : step : round(boundary_y(j)) + if ~BW_img(yy,boundary_x(j)) + flag = 1; + break + end + end + else + k = (boundary_y(i) - boundary_y(j))/(boundary_x(i) - boundary_x(j)); + step = sign(round(boundary_x(j))-round(boundary_x(i))); + for xx = round(boundary_x(i)) : step : round(boundary_x(j)) + yy = round(boundary_y(i) + k * (xx - boundary_x(i))); + if ~BW_img(yy,xx) + flag = 1; + break + end + end + end + if flag + continue + end + + chord_length = ((boundary_x(i) - boundary_x(j))^2 + (boundary_y(i) - boundary_y(j))^2).^0.5; + if chord_length > chord_length_max + chord_length_max = chord_length; + i_max = i; + j_max = j; + end + end +end +chord_max_x = [boundary_x(i_max) boundary_x(j_max)]; +chord_max_y = [boundary_y(i_max) boundary_y(j_max)]; + +if chord_max_x(1) == chord_max_x(2) + ang = 0; +elseif chord_max_y(1) == chord_max_y(2) + ang = pi/2; +else + ang = pi/2 - atan((chord_max_y(2) - chord_max_y(1))/(chord_max_x(2) - chord_max_x(1))); +end +BW_img_rotated = imrotate(BW_img,-ang/pi * 180,'crop'); +chord_max_x_R = round((chord_max_x - x0) * cos(ang) - (chord_max_y - y0) * sin(ang) + x0); +chord_max_y_R = round((chord_max_x - x0) * sin(ang) + (chord_max_y - y0) * cos(ang) + y0); +chord_max_x_R_mid = round(mean(chord_max_x_R)); +chord_max_ppd_max = -1; + +for yy = min(chord_max_y_R):max(chord_max_y_R) + + left = BW_img_rotated(yy,1:chord_max_x_R_mid); + x1 = find(left==0,1,'last'); + right = BW_img_rotated(yy,chord_max_x_R_mid:s(2)); + x2 = find(right==0,1,'first') + chord_max_x_R_mid; + if x2 - x1 > chord_max_ppd_max + chord_max_ppd_max = x2 - x1; + chord_max_ppd_max_x1 = x1; + chord_max_ppd_max_x2 = x2; + chord_max_ppd_max_yy = yy; + end +end + +chords_xy = cat(2,chord_max_x_R(:),chord_max_y_R(:),... + [chord_max_ppd_max_x2;chord_max_ppd_max_x1],[chord_max_ppd_max_yy;chord_max_ppd_max_yy]); +% catch +% keyboard +% end +end \ No newline at end of file diff --git a/cell_connect_tracks.m b/cell_connect_tracks.m new file mode 100644 index 0000000..35c46da --- /dev/null +++ b/cell_connect_tracks.m @@ -0,0 +1,134 @@ +function cell_connect_tracks(case_string) + +track_dir = fullfile('cases',case_string,'tracks'); + +num_frms = numel(dir([track_dir ,'/*.mat'])); + +tracks_euler = zeros(num_frms,1000,9); + +for volume = 1 : num_frms-1 + + load([track_dir,'/tracks_',num2str(volume,'%.4d'),'.mat']); + + tracks_euler(volume,1:size(tracks,1),:) = tracks; + +end; + +num_of_tracks = 0; + +tracks_lagrangian = zeros(num_frms,4,1000); + +for volume = 1 : num_frms + + for p = 1:size(tracks_euler,2) + + if tracks_euler(volume,p,1) ~= 0 + + num_of_tracks = num_of_tracks + 1; + + tracks_lagrangian(volume,:,num_of_tracks) = tracks_euler(volume,p,[1,4:6]); + tracks_lagrangian(volume+1,:,num_of_tracks) = tracks_euler(volume,p,[2,7:9]); + target_num = tracks_euler(volume,p,2); + tracks_euler(volume,p,:) = zeros(1,9); + + for nxt_v = volume + 1 : num_frms + + ind = tracks_euler(nxt_v,:,1) == target_num; + + if sum(ind) == 0 + break + end + + tracks_lagrangian(nxt_v,:,num_of_tracks) = tracks_euler(nxt_v,ind,[1,4:6]); + tracks_lagrangian(nxt_v+1,:,num_of_tracks) = tracks_euler(nxt_v,ind,[2,7:9]); + + target_num = tracks_euler(nxt_v,ind,2); + tracks_euler(nxt_v,ind,:) = zeros(1,9); + end + end + end + +end + +% display(sprintf('Number of tracks : %.4d',num_of_tracks)); + + +fst_frame = zeros(1, num_of_tracks); +lst_frame = zeros(1, num_of_tracks); + +for p = 1 : num_of_tracks + + current_track = tracks_lagrangian(:,:,p); + ind = find(current_track(:,1) ~=0); + fst_frame(p) = min(ind); + lst_frame(p) = max(ind); + +end; + + +for p1 = 1 : num_of_tracks + + p1_t = lst_frame(p1); + if p1_t == 0 + continue + end + + p1_p = tracks_lagrangian(p1_t ,2:4,p1); + + for p2 = 1 : num_of_tracks + + p2_t = fst_frame(p2); + if p2_t == 0 + continue + end + + p2_p = tracks_lagrangian(p2_t ,2:4,p2); + + dis = sum((p1_p-p2_p).^2)^0.5; + + if (p1_t == p2_t)&&(dis<5) + tracks_lagrangian(fst_frame(p2):lst_frame(p2),:,p1) = tracks_lagrangian(fst_frame(p2):lst_frame(p2) ,:,p2); + lst_frame(p1) = lst_frame(p2); + fst_frame(p2) = 0; + lst_frame(p2) = 0; + end; + + if (p1_t + 1 == p2_t)&&(dis<10) + tracks_lagrangian(fst_frame(p2):lst_frame(p2),:,p1) = tracks_lagrangian(fst_frame(p2):lst_frame(p2) ,:,p2); + tracks_lagrangian(lst_frame(p1)+1,:,p1) = ... + (tracks_lagrangian(lst_frame(p1),:,p1) + tracks_lagrangian(lst_frame(p1)+2,:,p1) )/2; + lst_frame(p1) = lst_frame(p2); + fst_frame(p2) = 0; + lst_frame(p2) = 0; + end; + + end +end + +ind = find(fst_frame ~=0); +tracks_lagrangian = tracks_lagrangian(:,:,ind); +num_of_tracks = length(ind); +% display(sprintf('Number of tracks after connection : %.4d',num_of_tracks)); + +duration = zeros(num_of_tracks,1); +position = zeros(num_of_tracks,3); +displacement = zeros(num_of_tracks,1); + +for p = 1 : num_of_tracks + + current_track = tracks_lagrangian(:,:,p); + + ind = find(current_track(:,1) ~= 0); + s0 = current_track(ind(1),2:4); + st = current_track(ind(end),2:4); + + duration(p) = length(ind); + position(p,:) = mean(current_track(ind,2:4),1); + displacement(p) = sum((st - s0).^2)^0.5; + +end + +% [duration,order] = sort(duration,'descend'); +[duration,order] = sort(displacement,'descend'); +tracks_order = tracks_lagrangian(:,:,order); +save(fullfile(track_dir,'tracks.mat'),'tracks_order') diff --git a/cell_core_segmentation.m b/cell_core_segmentation.m new file mode 100644 index 0000000..e907af9 --- /dev/null +++ b/cell_core_segmentation.m @@ -0,0 +1,223 @@ +function core = cell_core_segmentation(img,peak_min,con) + +% ========================================================================= +% segment cells via a 2D erosion\dilation method. +% by Tianqi Guo, Feb 2018 +% ========================================================================= + +% index-shift matrix for 4\8 pixel neighbourhood +if con == 4 + dx = [ 0 0 -1 1]; + dy = [ -1 1 0 0]; +elseif con == 8 + dx = [ 0 0 -1 1 -1 -1 1 1]; + dy = [-1 1 0 0 -1 1 -1 1]; +end + +% range for peak intensities +[s1,s2] = size(img); + +% range for peak intensities +peak_lo = peak_min; +peak_hi = max(img(:)); +peak_step = 5; + +% ========================================================================= +% erosion procedure + +img = double(img); +% count for discovered local intensity peaks and particle numbers +p = 0; + +level_set = peak_lo : peak_step : peak_hi + peak_step; +level_num = size(level_set,2); +region_peak_ind = cell(level_num,1); + +% gradually increase intensity threshold for peaks +for i = 1 : level_num + + level = level_set(i); + + % all intensities below current threshold are set to zero + img_l = img; + img_l(img_l < level) = 0; + + % find the connected regions for the current level + % each region potentially has a local maximum, i.e. a peak + CC = bwconncomp(img_l); + peak_ind = []; + for j = 1 : size(CC.PixelIdxList,2) + + region_ind = cell2mat(CC.PixelIdxList(j)); + [~,ind] = max(img(region_ind)); + peak_ind(j) = region_ind(ind(1)); + + end + + region_peak_ind{i} = peak_ind; +end + +% peaks with pixel index (ii, jj) +% peaks = zeros(particle_num,3); +peaks = []; + +for i = 2 : level_num + + peaks_old = cell2mat(region_peak_ind(i-1)); + peaks_new = cell2mat(region_peak_ind(i)); + + flagSubset = ~ismember(peaks_old,peaks_new); + + [ii, jj] = ind2sub(size(img),peaks_old(flagSubset)); + + dp = sum(flagSubset); + + peaks(p+1:p+dp,:) = [ii; jj]'; + + p = p + dp; + +end + +[ii, jj] = ind2sub(size(img),cell2mat(region_peak_ind(level_num))); +dp = length(ii); +peaks(p+1:p+dp,:) = [ii; jj]'; +p = p + dp; + +peaks = peaks(1:p,:); + +% ========================================================================= +% dilation procedure + +% the matrix that marks each voxel in particle cores with the particle number +core = zeros(size(img)); + +% a slower but more robust version of the dilation procedure +edge_i = zeros(1,p); +boundaries = cell(1,p); +for current_particle = 1 : p + core_i = peaks(current_particle,1); + core_j = peaks(current_particle,2); + boundaries{current_particle} = [core_i, core_j]; + edge_i(current_particle) = img(core_i, core_j) * 0.25; + core(core_i, core_j) = current_particle; +end + +while ~isempty(boundaries) + + new_boundaries = {}; + + for i = 1:length(boundaries) + + boundary = boundaries{i}; + new_boundary = []; + + for j = 1:size(boundary,1) + + x = boundary(j,1); + y = boundary(j,2); + current_particle = core(x,y); + + for search_voxel = 1 : size(dx,2) + + % the indices of the surrounding voxel + xx = x + dx(search_voxel); + yy = y + dy(search_voxel); + + if (xx < 1) || (xx > s1) || (yy < 1) || (yy > s2) || (core(xx,yy) > 0) + continue + end; + + + % if this voxel is brighter than the basic threshold + % and dimmer than the current voxel + if (img(xx,yy) >= edge_i(current_particle)) && (img(xx,yy) <= img(x,y)) %img(core_i,core_j)) % + + % mark this voxel with current particle number + core(xx,yy) = current_particle; + new_boundary = [new_boundary;xx yy]; + + end + end + + + end + if ~isempty(new_boundary) + new_boundaries{length(new_boundaries)+1} = new_boundary; + end + end + boundaries = new_boundaries; + +end + + +% a faster but less robust version of the dilation procedure + +% % for each peak recorded in the erosion procedure +% for current_particle = 1 : p +% +% % initialize the flood-fill with the peak location as the seed +% seed = peaks(current_particle,:); +% core_i = seed(1); +% core_j = seed(2); +% core(core_i,core_j) = current_particle; +% edge_i = img(core_i,core_j) * 0.50; +% +% % place the seed into the search stack +% l = 1; +% list = zeros(200,2); +% list(l,:) = seed; +% +% % the flood-fill procedure that marks the voxels with particle number +% while l > 0 +% +% % take the voxel from the tail of the stack +% pixel = list(l,:); +% x = pixel(1); +% y = pixel(2); +% l = l - 1; +% +% % look around the current voxel in the neighbors +% for search_voxel = 1 : size(dx,2) +% +% % the indices of the surrounding voxel +% xx = x + dx(search_voxel); +% yy = y + dy(search_voxel); +% +% if (xx < 1) || (xx > s1) || (yy < 1) || (yy > s2) +% continue +% end; +% +% % if this voxel is brighter than the basic threshold +% % and dimmer than the current voxel +% if (img(xx,yy) >= edge_i) && (img(xx,yy) <= img(x,y)) % img(core_i,core_j)) % +% +% % if this voxel has not been visited before +% if (core(xx,yy) == 0) +% +% % mark this voxel with current particle number +% core(xx,yy) = current_particle; +% +% % place this voxel into the stack +% l = l + 1; +% list(l,:) = [xx,yy]; +% +% end +% +% % if this voxel has been visited from other cores +% if (core(xx,yy) > 0) && (core(xx,yy) ~= current_particle) +% +% % this voxel doesn't belong to a particle core +% % since it is under the influece of more than one particle +% core(xx,yy) = -1; +% +% end +% end +% end +% end +% +% end +% +% % mark all the non-core voxels with zero +% core(core == -1) = 0; +% +% end \ No newline at end of file diff --git a/cell_downsample_Hessian_segmentation.m b/cell_downsample_Hessian_segmentation.m new file mode 100644 index 0000000..54c6132 --- /dev/null +++ b/cell_downsample_Hessian_segmentation.m @@ -0,0 +1,158 @@ +function cell_downsample_Hessian_segmentation(props) + +% ========================================================================= +img_fname = props.img_fname; +case_string = props.case_string; +optionFigure = props.optionFigure; +window_range = props.window_range; +downsample_factor = props.downsample_factor; +con = props.connectivity; + +peak_min = props.peak_min; +basic_t = props.basic_t; +area_min = props.area_min; +area_max = props.area_max; +% ========================================================================= + +info = imfinfo(img_fname); +num_images = numel(info); + +x1 = window_range(1); +x2 = window_range(2); +y1 = window_range(3); +y2 = window_range(4); + +downsize_halfWin = max(floor(downsample_factor/2),1); +ave_filter = fspecial('average', [downsample_factor downsample_factor]); + +case_dir = fullfile('cases',case_string); +mkdir(case_dir); +rmdir(case_dir,'s'); +mkdir(case_dir); + +out_dir_images = fullfile(case_dir,'images'); +mkdir(out_dir_images); + +out_dir_segmentation = fullfile(case_dir,'segmentation'); +mkdir(out_dir_segmentation); + +% if optionFigure +% +% out_dir_frames = fullfile(case_dir,'frames'); +% mkdir(out_dir_frames); +% +% myVideo = VideoWriter(fullfile(out_dir_frames,sprintf('undersampled_%.2d.avi',downsample_factor))); +% myVideo.FrameRate = 5; % Default 30 +% myVideo.Quality = 75; % Default 75 +% open(myVideo) +% figure('color','w','position',[0 200 1800 1000]) +% end +% + +parfor current_frame = 1:num_images + + img_raw = imread(img_fname, current_frame); + img_cropped = medfilt2(img_raw(y1:y2,x1:x2)); + + Y = prctile(img_cropped(:),99.9); + img_adjusted = uint8(double(img_cropped)/double(Y) * 0.9 * double(intmax('uint8'))); + + img_aveFiltered = uint8(filter2(ave_filter, img_adjusted)); + s = size(img_aveFiltered); + + n_x = floor(s(2)/downsample_factor); + n_y = floor(s(1)/downsample_factor); + img_aveFiltered = img_aveFiltered(1:n_y*downsample_factor,1:n_x*downsample_factor); + img_downSampled = img_aveFiltered(downsize_halfWin:downsample_factor:end-downsize_halfWin,downsize_halfWin:downsample_factor:end-downsize_halfWin); + + [lambda1,lambda2] = Hessian_2D(img_downSampled); + He_mask = (lambda2<0)&(lambda1<0); + He_mask([1:3,end-3:end],:) = 0; + He_mask(:,[1:3,end-3:end]) = 0; + cen_He = cell2mat(struct2cell(regionprops(He_mask,'Centroid'))'); + N_mask = length(cen_He); + + img_mean = ceil(mean(img_downSampled(:))); + BW = imbinarize(img_downSampled, img_mean/255); + cen_BW = cell2mat(struct2cell(regionprops(BW,'Centroid'))'); + [L, N_BW] = bwlabel(BW); + for r = 1:N_mask + [dis_min, true_num] = finding_true_num_2D(cen_BW,cen_He(r,:)); + L(L==true_num) = N_BW + 1; + end + BW_new = (L == N_BW + 1); + L_new = bwlabel(BW_new,8); + + img_HessianFiltered = img_downSampled; + img_HessianFiltered(~BW) = 0; + + core = cell_core_segmentation(img_HessianFiltered,peak_min,con); + + seg_raw = repelem(core,downsample_factor,downsample_factor); + seg = cell_refinement(seg_raw,img_aveFiltered,basic_t,area_min,area_max); + parsave(fullfile(out_dir_segmentation,sprintf('segmentation_%.4d',current_frame)),seg) + + img_adjusted = img_adjusted(1:n_y*downsample_factor,1:n_x*downsample_factor); + filename = fullfile(out_dir_images,sprintf('image_%.4d.tif',current_frame)); + imwrite(img_adjusted,filename); %for saving preprocessed images + +% if optionFigure +% +% subplot(361) +% imshow(img_cropped) +% title('Raw image') +% +% subplot(362) +% imshow(img_adjusted) +% title('Intensity adjusted') +% +% subplot(363) +% imshow(img_aveFiltered) +% title('Moving averaging filtered') +% +% subplot(367) +% imshow(img_downSampled) +% title(sprintf('Downsampled by factor of %.1d',downsample_factor)) +% +% subplot(368) +% imagesc(lambda1) +% caxis([-10 10]) +% axis equal off +% title('Hessian \lambda_1') +% +% subplot(369) +% imagesc(lambda2) +% axis equal off +% caxis([-10 10]) +% title('Hessian \lambda_2') +% +% subplot(3,6,13) +% imshow(img_HessianFiltered) +% title('Hessian filtered') +% +% subplot(3,6,14) +% imagesc(core) +% axis equal off +% title('Segmented cell core') +% +% subplot(3,6,15) +% imagesc(seg) +% axis equal off +% title('Refined cell profile') +% +% subplot(3,6,[4 5 6 10 11 12 16 17 18]) +% imshow(img_adjusted); hold on; +% contour(seg~=0,1,'b','linewidth',2); hold on; +% % plot(cell_pos(:,1),cell_pos(:,2),'rX','markersize',5,'linewidth',2,'linesmoothing','on'); hold on; +% % ellipse(cell_paras(:,1),cell_paras(:,2),cell_paras(:,3),cell_pos(:,1),cell_pos(:,2),'r',50); hold on; +% title(sprintf('Segmented cell from frame = %.3d',current_frame)) +% +% % print(fullfile(out_dir_frames,sprintf('segmentation_frame=%.4d.tiff',current_frame)),'-dtiff','-r100') +% set(gca,'nextplot','replacechildren'); +% frame = getframe(gcf); +% writeVideo(myVideo, frame) +% end +end +% if optionFigure +% close(myVideo) +% end diff --git a/cell_parameterization.m b/cell_parameterization.m new file mode 100644 index 0000000..a609245 --- /dev/null +++ b/cell_parameterization.m @@ -0,0 +1,166 @@ +function cell_parameterization(props) + +case_string = props.case_string; +optionDebug = props.optionDebug; + +case_dir = fullfile('cases',case_string); + +if optionDebug + debug_dir = fullfile(case_dir,'debug'); + mkdir(debug_dir) +end + +seg_dir = fullfile(case_dir,'segmentation'); +seg_list = dir(fullfile(seg_dir,'*.mat')); + +GSM_dir = fullfile(case_dir,'GSM'); +mkdir(GSM_dir); + +particle_dir = fullfile(case_dir,'particles'); +mkdir(particle_dir); + +image_dir = fullfile(case_dir,'images'); +img_list = dir(fullfile(image_dir,'*.tif')); +num_images = length(img_list); + +w = 100; + +measures = cell(num_images,1); + +for current_frame = 1:num_images + + img_raw = imread(fullfile(image_dir,img_list(current_frame).name)); + + load(fullfile(seg_dir,seg_list(current_frame).name)); + seg = data; + + seg_pad = padarray(seg,[w w],'both'); + img_pad = padarray(img_raw,[w w],'both'); + + cell_N = max(seg(:)); + measures_currentFrame = zeros(cell_N,27); + + parfor j = 1:cell_N + + BW_j = (seg_pad == j); + + img_j = img_pad; + img_j(~BW_j) = 0; + + cen = round(cell2mat(struct2cell(regionprops(BW_j,'Centroid')))); + + if isempty(cen) + continue + end + + [y,x] = ind2sub(size(BW_j),find(BW_j)); + ww = max(max(abs(y - cen(2))), max(abs(x - cen(1)))) + 25; + + x1 = cen(1)-ww; + x2 = cen(1)+ww; + y1 = cen(2)-ww; + y2 = cen(2)+ww; + + img_crop = double(img_j(y1:y2,x1:x2)); + BW_crop_edge = edge(BW_j(y1:y2,x1:x2),'log'); + BW_crop = imfill(BW_crop_edge,'holes'); + area = cell2mat(struct2cell(regionprops(BW_crop,'area'))); + if length(area) > 1 + [~,max_ind] = max(area); + BW_crop = (bwlabel(BW_crop) == max_ind); + end + + % figure('position',[0 0 900 900]) + % subplot(231) + % imshow(img_raw_crop) + % subplot(234) + % imagesc(img_crop); axis equal off + % subplot(232) + % imshow(BW_crop) + % subplot(233) + % imshow(BW_crop_edge) + % subplot(232) + % imshow(BW_crop) + % subplot(233) + % imshow(BW_crop_edge) + + measure = calculate_GSM(img_crop,BW_crop,optionDebug); + measure(1) = measure(1) + x1 - w; + measure(2) = measure(2) + y1 - w; + measure(3) = measure(3) + x1 - w; + measure(4) = measure(4) + y1 - w; + if find(inf(isnan(measure) | isinf(measure))) + continue + end + measures_currentFrame(j,:) = measure; + +% if optionDebug +% +% [y,x] = ind2sub(size(BW_crop_edge),find(BW_crop_edge)); +% +% subplot(241) +% imshow(img_crop); hold on; +% plot(x,y,'.r'); +% title('MFI image') +% +% +% print(fullfile(debug_dir,sprintf('_%.02d.tif',j)),'-dtiff','-r200') +% end + + end + ind = measures_currentFrame(:,1) == 0; + measures_currentFrame(ind,:) = []; + save(fullfile(GSM_dir,sprintf('%0.4d',current_frame)),'measures_currentFrame') + measures{current_frame} = measures_currentFrame; +end + +save(fullfile(case_dir,'GSM_all'),'measures') +measures_all = cell2mat(measures); + +cats = {'IWC','Center of mass',... + 'I_{mean}','I_{max}','Area (pixel^2)', 'Perimeter (pixel)',... + 'Major axis (pixel)', 'Minor axis (pixel)','Equiv R (pixel)', 'Orientation (°)', ... + 'Elongation (FFT)','Ecc (chord)','Ecc (ellipse)', 'Elongatedeness', 'Ratio (Rb/Ra)',... + 'Circ (AP)', 'Circ (r_{std})', 'Circ (r_{avg})', 'Circ (r_{max})', ... + 'Ellipticity','Squareness', 'Triangularity', 'Rectangularity', 'Convexity (A)', 'Convexity (P)'}; +sub_pos = [0 0 0 0 1 2 3 4 7 8 9 10 13 14 15 16 17 19 20 21 22 25 26 27 28 29 30]; +close all; +% figure('color','white','position',[0 0 1800 1000]) +% for i = 5:27 +% subplot(5,6,sub_pos(i)) +% measure = measures_all(:,i); +% ind = find(isnan(measure) | isinf(measure)); +% measure(ind) = []; +% measure_avg = abs(mean(measure)); +% measure_std = var(measure).^0.5; +% edges = linspace(measure_avg-3*measure_std,measure_avg+3*measure_std,20); +% histogram(measure,edges,'Normalization','probability') +% xlim([min(measure) max(measure)]) +% title({cats{i-2},sprintf('RMS/MEAN = %.2f',measure_std/measure_avg)}) +% xlim([measure_avg-3*measure_std measure_avg+3*measure_std]) +% end +% print(fullfile(case_dir,'GSM_pdf.tif'),'-dtiff','-r300') + + +measures_ind = [10,9,12,6,7,16,19,23]; +for current_frame = 1:num_images + + measures_currentFrame = measures{current_frame}; + N = length(measures_currentFrame); + pos = cat(2,measures_currentFrame(:,1:2),zeros(size(measures_currentFrame,1),1)); + paras = measures_currentFrame(:,measures_ind); + paras(:,3) = pi/2 - paras(:,3)/180 * pi; + particles = []; + particles.pos = pos; + particles.paras = paras; + save(fullfile(particle_dir,sprintf('%0.4d',current_frame)),'particles') + +end + +% +% +% particles = []; +% particles.pos = cell_pos; +% particles.paras = cell_paras; +% file_name = fullfile(out_dir_particles,sprintf('particles_%.4d.mat',current_frame)); +% parsave(file_name,particles) diff --git a/cell_plot_sequence.m b/cell_plot_sequence.m new file mode 100644 index 0000000..4d6edd8 --- /dev/null +++ b/cell_plot_sequence.m @@ -0,0 +1,78 @@ +function cell_plot_sequence(case_string,time_step) + +case_dir = fullfile('cases',case_string); +% time_step = 1; +% case_string = 'ca1a_n2'; +% case_dir = '/home/shannon/b/aether/Projects/Tumor_Microenvironment/Analysis/Brian/Spring_2019/03032019/cases/ca1hFN3_n1_bsubtracted_v1'; + +img_dir = fullfile(case_dir,'images'); +track_dir = fullfile(case_dir,'tracks'); +particle_dir = fullfile(case_dir,'particles'); + +image_list = dir(fullfile(img_dir,'*.tif')); +num_imgs = length(image_list); + +particle_list = dir(fullfile(particle_dir,'*.mat')); + +load(fullfile(track_dir,'tracks.mat')) + +out_dir = 'results'; +mkdir(out_dir); + +myVideo = VideoWriter(fullfile(out_dir,[case_string,'.avi'])); +myVideo.FrameRate = 30; % Default 30 +myVideo.Quality = 50; % Default 75 +open(myVideo) + +num_tracks = size(tracks_order,3); +% colors = distinguishable_colors(num_tracks,'k'); +colors = rand(num_tracks,3); +figure('position',[0 0 900 900],'InvertHardCopy', 'off','color','w') + +for i = 1: num_imgs - 2 + time = (i - 1) * time_step; + file_name = fullfile(img_dir,image_list(i).name); + img_raw = imread(file_name); + img = medfilt2(img_raw,[3 3]); + + current_track = reshape(tracks_order(i,:,:),4,num_tracks)'; + ind = find(current_track(:,1)~=0); + + %[rb(ind) ra(ind) ang(ind) area(ind) I(ind)]; + + file_name = fullfile(particle_dir,particle_list(i).name); + load(file_name) + + pos = particles.pos; + paras = particles.paras; + + imshow( img_raw,'initialmagnification','fit'); hold on; + + for j = 1:length(ind) + track_ID = ind(j); + track_j = tracks_order(:,:,track_ID); + ind_t = find(track_j(:,1) ~=0); + times = ind_t(1):i+1; + x_t = track_j(times,2); + y_t = track_j(times,3); + ind_t = (x_t.*y_t ~= 0); + x_t = x_t(ind_t); + y_t = y_t(ind_t); + particle_id = current_track(track_ID,1); + plot(x_t(1),y_t(1),'color',colors(track_ID,:),'marker','X','markersize',5,'linewidth',2,'linesmoothing','on'); hold on; + %ellipse(paras(particle_id,1),paras(particle_id,2),paras(particle_id,3),pos(particle_id,1),pos(particle_id,2),colors(track_ID,:),50); hold on; + plot(x_t,y_t,'-','color',colors(track_ID,:),'linewidth',2,'linesmoothing','on'); hold on; + + end + hold off; + title(sprintf('Time = %.2f frame',time)) + % file_name = [out_dir,sprintf('/t=%.4d.tif',i)]; + % print(file_name,'-dtiff','-r100') + % frame = imread(file_name); + + set(gca,'nextplot','replacechildren'); + frame = getframe(gcf); + writeVideo(myVideo, frame) + +end +close(myVideo); \ No newline at end of file diff --git a/cell_postprocessing.m b/cell_postprocessing.m new file mode 100644 index 0000000..9502fe5 --- /dev/null +++ b/cell_postprocessing.m @@ -0,0 +1,271 @@ +function cell_postprocessing(props) + +plot_dXdY_max = props.plot_dXdY_max; +plot_UV_max = props.plot_UV_max; +plot_cell_num_max = props.plot_cell_num_max; +d_acc_threshold = props.d_acc_threshold; +time_step = props.time_step; +case_string = props.case_string; + +track_dir = fullfile('cases',case_string,'tracks'); +particle_dir = fullfile('cases',case_string,'particles'); +plot_dir = 'results'; +mkdir(plot_dir) +load([track_dir,'/tracks.mat']) +num_frms = numel(dir(fullfile(particle_dir ,'*.mat'))); +num_tracks = size(tracks_order,3); +close all +figure('position',[0 0 1000 1000],'color','w') + +dX = []; +dY = []; +U = []; +V = []; +d_euclid = zeros(1,num_tracks); +d_accum = zeros(1,num_tracks); +xy_end = zeros(num_tracks,2); +c = 0; +flag = zeros(1,num_tracks); +motile_num = zeros(1,num_frms); +longest_length = 0; + +for p = 1 : num_tracks + + % extract current Lagrangian track + current_track = tracks_order(:,:,p); + ind = current_track(:,2) ~= 0; + longest_length = max(longest_length,sum(ind)); + current_track = current_track(ind,:); + + % subtract the original place + X = current_track(:,2) - current_track(1,2); + Y = current_track(:,3) - current_track(1,3); + + % instantaneous displacement over each frame + U_p = X(2:end) - X(1:end-1); + V_p = Y(2:end) - Y(1:end-1); + d_accum_p = sum((U_p.^2 + V_p.^2).^0.5); + + if d_accum_p >= d_acc_threshold + frames = find(ind); + try + for f = frames + motile_num(f) = motile_num(f) + 1; + end + catch + keyboard + end + c = c + 1; + flag(p) = 1; + d_accum(c) = d_accum_p; + U = [U; U_p(:)]; + V = [V; V_p(:)]; + + % integral displacement from the original place + dX_p = X(end); + dY_p = Y(end); + xy_end(c,:) = [dX_p dY_p]; + d_euclid(c) = (dX_p^2 + dY_p^2)^0.5; + dX = [dX; dX_p(:)]; + dY = [dY; dY_p(:)]; + end +end + +tracks_thresholded = tracks_order(:,:,find(flag)); +save([track_dir,'/tracks_thresholded.mat'],'tracks_thresholded') + +ind = 1:c; +d_euclid = d_euclid(ind); +d_accum = d_accum(ind); +xy_end = xy_end(ind,:); +% Directness +D_p = d_euclid./d_accum; +D_a = mean(D_p); + +% Center of Mass +M_end = mean(xy_end,1); + +% Forward Migration Index +y_FMI = mean(xy_end(:,2)./d_accum(:)); +x_FMI = mean(xy_end(:,1)./d_accum(:)); + +time_axis = (0:num_frms-1) * time_step/60; + +cell_num = zeros(num_frms,1); +file_list = dir(fullfile(particle_dir,'*.mat')); +for t = 1: num_frms + file_name = fullfile(particle_dir,file_list(t).name); + load(file_name); + cell_num(t) = length(particles.pos); +end + +results = []; +results.time_step = time_step; +results.case_string = case_string; +results.d_acc_threshold = d_acc_threshold; + +% results.mean_int_threshold = basic_t; +% results.min_peak_int_threshold = peak_min; +% results.min_area = area_min; + +results.U = U; +results.V = V; +results.dX = dX; +results.dY = dY; +results.cell_num = cell_num; +results.motile_num = motile_num; + +results.d_accum = d_accum; % added by Brian +results.D_a = D_a; +results.M_end = M_end; +results.FMI = [x_FMI,y_FMI]; + +save(fullfile(plot_dir,case_string),'results') +%========================================================================== +subplot(339) +plot(0,0) +text(0,8,['Case: ',strrep(case_string,'_',' ')]) +text(0,6,sprintf('Averaged directness = %.2f',D_a)) +text(0,4,sprintf('Averaged center of Mass = (%.2f, %.2f)',M_end(1),M_end(2))) +text(0,2,sprintf('Forward Migration Index = (%.2f, %.2f)',x_FMI,y_FMI)) +axis equal off +xlim([0 10]) +ylim([0 10]) + +subplot(338) +cell_num_s = smooth(cell_num,0.2); +motile_num_s = smooth(motile_num,0.2); +plot(time_axis,cell_num_s,'O','linewidth',2,'markersize',5); hold on; +plot(time_axis,motile_num_s,'O','linewidth',2,'markersize',5); hold on; +xlim([0 time_axis(end)]) +ylim([0 plot_cell_num_max]) +title('Cell counts over time') +grid on; +xlabel('t (hours)') +ylabel('Cell counts') +legend('Total cell number','Motile cell number','location','southoutside') + +%========================================================================== + + +% histogram +% dXdY_max = round(max(abs([dX(:);dY(:)])),-1); +% bin_size = 10; +% bins_half = (0:bin_size:dXdY_max) + bin_size/2; +% bins = [bins_half - (dXdY_max + bin_size),bins_half]; +% [xx,yy] = meshgrid(edges,edges); +edges = - plot_dXdY_max :10: plot_dXdY_max; +N = hist3([dX,dY],'Ctrs',{edges,edges})'; +N_x = sum(N,1); +N_y = sum(N,2); + +subplot(331) +bar(edges,N_x/sum(N_x),'EdgeColor','none') +xlim([-1 1]*plot_dXdY_max) +title('Histogram of (X_f-X_0)') +xlabel('X_f-X_0 (pixel)') +grid on; + +subplot(333) +bar(edges,N_y/sum(N_y),'EdgeColor','none') +xlim([-1 1]*plot_dXdY_max) +title('Histogram of (Y_f-Y_0)') +xlabel('Y_f-Y_0 (pixel)') +grid on; + +subplot(332) +% color by tracks +colors_1 = distinguishable_colors(c); + +% color by time - specific colormap +colors_2 = flipud(jet(longest_length)); + +% color by time - manual colormap +cc = (longest_length:-1:1)/longest_length; +colors_3 = [cc(:),cc(:),ones(longest_length,1)]; + +colors = colors_2; + +p_set = find(flag); +for p = 1:20 %c + + % extract current Lagrangian track + current_track = tracks_order(:,:,p_set(p)); + ind = current_track(:,2) ~= 0; + current_track = current_track(ind,:); + + % subtract the original place + X = current_track(:,2) - current_track(1,2); + Y = current_track(:,3) - current_track(1,3); + + % color by time - specific colormap + for t = 2:length(X) + plot([X(t) X(t-1)],[Y(t) Y(t-1)],'color',[colors(t,:),0.5],'linewidth',2,'linesmoothing','on'); hold on; + end + plot(X(end),Y(end),'Ob',... + 'markerfacecolor','b','linewidth',2,'linesmoothing','on','markersize',2); hold on; + +% color by tracks +% plot(X,Y,'color',colors(p,:),'linewidth',1,'linesmoothing','on'); hold on; +% plot(X(end),Y(end),'O','color',colors(p,:),... +% 'markerfacecolor',colors(p,:),'linewidth',2,'linesmoothing','on','markersize',2); hold on; +end +plot([0 0],[-1 1]*plot_dXdY_max,'--r','linewidth',2); hold on; +plot([-1 1]*plot_dXdY_max,[0 0],'--r','linewidth',2); +xlabel('X - X_0 (pixel)') +ylabel('Y - Y_0 (pixel)') +axis equal +xlim([-1 1]*plot_dXdY_max) +ylim([-1 1]*plot_dXdY_max) +title('Cell movement from original place') + + +% UV_max = round(max(abs([U(:);V(:)])),-1); +% bin_size = 1; +% bins_half = (0:bin_size:UV_max) + bin_size/2; +% bins = [bins_half - (UV_max + bin_size),bins_half]; +% [xx,yy] = meshgrid(bins,bins); +edges = - plot_UV_max : plot_UV_max; +N = hist3([U,V],'Ctrs',{edges,edges})'; +N_x = sum(N,1); +N_y = sum(N,2); + +n = props.roseBins; +[theta,rho] = cart2pol(U,V); +r = (0:n)'/n * plot_UV_max * sqrt(2); +t = pi*(-n:n)/n; +X = r*cos(t); +Y = r*sin(t); +N = hist3([theta,rho],'edges',{t,r+max(rho)/n})'; + +subplot(335) +% h = pcolor(xx-bin_size/2,yy-bin_size/2,N); hold on +h = pcolor(X,Y,N); colorbar; axis equal tight; hold on +%set(h, 'EdgeColor', 'none') +%shading interp +plot([0 0],[-200 200],'--r','linewidth',2); hold on; +plot([-200 200],[0 0],'--r','linewidth',2); +xlabel('U (pixel/frame)') +ylabel('V (pixel/frame)') +axis equal +xlim([-1 1]*plot_UV_max) +ylim([-1 1]*plot_UV_max) +title({'2D histogram',' of movement per frame'}) +hold on +colorbar + +subplot(334) +bar(edges,N_x/sum(N_x),'EdgeColor','none') +xlim([-1 1]*plot_UV_max) +title('Histogram of U') +xlabel('U (pixel/frame)') +grid on; + +subplot(336) +bar(edges,N_y/sum(N_y),'EdgeColor','none') +xlim([-1 1]*plot_UV_max) +title('Histogram of V') +xlabel('V (pixel/frame)') +grid on; +print(fullfile(plot_dir,case_string),'-dtiff','-r300') +%========================================================================== \ No newline at end of file diff --git a/cell_processing_main.m b/cell_processing_main.m new file mode 100644 index 0000000..5d3eb02 --- /dev/null +++ b/cell_processing_main.m @@ -0,0 +1,136 @@ +clc; +clear; +close all; +warning('off','all') +set(0,'defaultAxesFontSize',10) +set(0,'defaultAxesFontName', 'DejaVu Sans') +set(0,'defaultTextFontName', 'DejaVu Sans') + +addpath(genpath('common_codes')); +% ========================================================================= + +raw_data_dir = uigetdir('ADD YOUR FOLDER DIRECTORY CONTAIING THE STACKED TIFF FILES'); +disp(['Data directory: ',raw_data_dir]) + +case_list = dir(fullfile(raw_data_dir,'*.tif')); + +N_cases = length(case_list); + +case_set = 1:N_cases; + +optionProcessing = 1; +optionPostprocessing = 0; +optionVideo = 1; + +if optionProcessing + plot_dir = 'results'; + mkdir(plot_dir) + window_range = zeros(N_cases,4); + + basic_t = zeros(N_cases,1); + peak_min = zeros(N_cases,1); + area_min = zeros(N_cases,1); + area_max = zeros(N_cases,1); + time_step = zeros(N_cases,1); + d_acc_threshold = zeros(N_cases,1); + + + for i = case_set + close all + figure('position',[200 100 800 800]) + C = strsplit(case_list(i).name,'.'); + case_string = strrep(C{1},'_',' '); + img_fname = fullfile(case_list(i).folder,case_list(i).name);%fullfile(raw_data_dir,case_list(i).name); + img = imread(img_fname, 1); + imagesc(img); hold on; + axis equal + title(['Select region for case: ',case_string]); + [x,y] = ginput(2); + x0 = min(round(x)); + x1 = max(round(x)); + y0 = min(round(y)); + y1 = max(round(y)); + rectangle('Position',[x0 y0 x1-x0 y1-y0],'EdgeColor','r','LineWidth',3,'linestyle','--') + print(fullfile(plot_dir,[case_string,'_ROI']),'-dtiff','-r300') + window_range(i,:) = [x0 x1 y0 y1]; + + close all + prompt = {'Basic intensity thresholding:',... + 'Minimum peak intensity:',... + 'Minimum cell area:',... + 'Maximum cell area:',... + 'Time step (min) between frames:',... + 'd_acc to threshold motile cells:'}; + tt = 'Processing parameters'; + dims = [1 35]; + definput = {'20','25','3','1000','30','50'}; + answer = inputdlg(prompt,tt,dims,definput); + basic_t(i) = str2double(answer{1}); + peak_min(i) = str2double(answer{2}); + area_min(i) = str2double(answer{3}); + area_max(i) = str2double(answer{4}); + time_step(i) = str2double(answer{5}); + d_acc_threshold(i) = str2double(answer{6}); + + end + + + props = []; + props.optionFigure = 0; + props.downsample_factor = 7; + + for i = case_set + + C = strsplit(case_list(i).name,'.'); + case_string = C{1}; + disp(['Processing case: ',case_string]) + + props.case_string = case_string; + props.window_range = window_range(i,:); + props.img_fname = fullfile(case_list(i).folder,case_list(i).name); + props.peak_min = peak_min(i); + props.basic_t = basic_t(i); + props.area_min = area_min(i); + props.area_max = area_max(i); + props.optionDebug = 0; + props.connectivity = 4; + + cell_downsample_Hessian_segmentation(props); + cell_parameterization(props); + cell_tracking(case_string); + cell_connect_tracks(case_string); + + end +end + +if optionPostprocessing + + props = []; + props.roseBins = 8; + props.plot_UV_max = 10; + props.plot_dXdY_max = 200; + props.plot_cell_num_max = 1000; + + for i = case_set + close all + C = strsplit(case_list(i).name,'.'); + props.case_string = C{1}; + + props.time_step = time_step(i); + props.d_acc_threshold = d_acc_threshold(i); + cell_postprocessing(props); + + end +end + +if optionVideo + time_step = 1; + for i = case_set + close all + C = strsplit(case_list(i).name,'.'); + case_string = C{1}; + + cell_plot_sequence(case_string,time_step); + + end +end \ No newline at end of file diff --git a/cell_refinement.m b/cell_refinement.m new file mode 100644 index 0000000..00eec77 --- /dev/null +++ b/cell_refinement.m @@ -0,0 +1,129 @@ +function seg = cell_refinement(seg_raw,img,basic_t,area_min,area_max) + +% dx = [0 0 0 -1 1 -1 -1 1 1]; +% dy = [0 -1 1 0 0 -1 1 1 -1]; + +dx = [0 0 0 -1 1 ]; +dy = [0 -1 1 0 0 ]; + +img = double(img); +s = size(seg_raw); +cell_N = max(seg_raw(:)); + +I_peaks = zeros(1,cell_N); +p_peaks = zeros(cell_N,2); + +for current_cell = 1:cell_N + + pixels = find(seg_raw == current_cell); + if isempty(pixels) + continue + end + [I_m,I_m_ind] = max(img(pixels)); + [I_m_r,I_m_c] = ind2sub(s,pixels(I_m_ind)); + + I_peaks(current_cell) = I_m; + p_peaks(current_cell,:) = [I_m_r,I_m_c]; + +end + +[~, ord] = sort(I_peaks,'ascend'); + +seg = zeros(size(img)); + +for p = 1:cell_N + + current_cell = ord(p); + + if I_peaks(current_cell) == 0 + continue + end + + edge_I = max(I_peaks(current_cell) * 0.5,basic_t); + + list = zeros(1000,2); + + pixels = find(seg_raw == current_cell); +% seg(pixels) = current_cell; + [rr, cc] = ind2sub(s,pixels); + + l = length(rr); + list(1:l,:) = horzcat(rr(:),cc(:)); + + + % the flood-fill procedure that marks the voxels with particle number + while l > 0 + + % take the voxel from the tail of the stack + pixel = list(l,:); + x = pixel(1); + y = pixel(2); + l = l - 1; + + % look around the current voxel in the neighbors + for search_voxel = 1 : size(dx,2) + + % the indices of the surrounding voxel + xx = x + dx(search_voxel); + yy = y + dy(search_voxel); + + if (xx < 1) || (xx > s(1)) || (yy < 1) || (yy > s(2)) + continue + end; + + % if this voxel is brighter than the basic threshold + % and dimmer than the current voxel + if (img(xx,yy) >= edge_I) && (img(xx,yy) <= img(x,y)) %I_peaks(current_cell)) %<= img(x,y)) % + + % if this voxel has not been visited before + if (seg(xx,yy) == 0) + + % mark this voxel with current particle number + seg(xx,yy) = current_cell; + + % place this voxel into the stack + l = l + 1; + list(l,:) = [xx,yy]; + + end + + % if this voxel has been visited from other cores + if (seg(xx,yy) > 0) && (seg(xx,yy) ~= current_cell) + + % this voxel doesn't belong to a particle core + % since it is under the influece of more than one particle + seg(xx,yy) = -1; + + end + end + end + end + +end + +% mark all the non-core voxels with zero +seg(seg == -1) = 0; +for current_cell = 1 : cell_N + + BW = imfill((seg == current_cell),'holes'); + seg(find(BW)) = 0; + + L = bwlabel(BW); + areas = cell2mat(struct2cell(regionprops(BW,'area'))); + [~,max_ind] = max(areas); + + if isempty(max_ind) + continue + end + BW = (L == max_ind(1)); + + ind = find(BW); + pixels_count = length(ind); + + if (pixels_count >= area_min) && (pixels_count <= area_max) + seg(find(BW)) = current_cell; + end + +end + + diff --git a/cell_tracking.m b/cell_tracking.m new file mode 100644 index 0000000..6d9af72 --- /dev/null +++ b/cell_tracking.m @@ -0,0 +1,32 @@ +function cell_tracking_v2(case_string) + +input_dir = fullfile('cases',case_string,'particles'); +frame_list = dir(fullfile(input_dir,'*.mat')); +n_frame = length(frame_list); + +s_radius = 10; +weights = [1,1,1,1,1,1,1,1,1]; + +out_dir = fullfile('cases',case_string,'tracks'); +mkdir(out_dir); +rmdir(out_dir,'s'); +mkdir(out_dir); + +for i = 1 : n_frame-1 + + file_name = fullfile(input_dir,frame_list(i).name); + load(file_name); + p1_pos = particles.pos; + p2_est = p1_pos; + paras_1 = particles.paras; + + file_name = fullfile(input_dir,frame_list(i+1).name); + load(file_name); + p2_pos = particles.pos; + paras_2 = particles.paras; + + tracks = weighted_nearest_neighbor3D_MP(p2_est,p1_pos,p2_pos,paras_1,paras_2,s_radius,weights); + file_name = fullfile(out_dir,sprintf('tracks_%.4d.mat',i)); + save(file_name,'tracks'); + +end \ No newline at end of file diff --git a/distinguishable_colors.m b/distinguishable_colors.m new file mode 100644 index 0000000..f9e9e25 --- /dev/null +++ b/distinguishable_colors.m @@ -0,0 +1,152 @@ +function colors = distinguishable_colors(n_colors,bg,func) +% DISTINGUISHABLE_COLORS: pick colors that are maximally perceptually distinct +% +% When plotting a set of lines, you may want to distinguish them by color. +% By default, Matlab chooses a small set of colors and cycles among them, +% and so if you have more than a few lines there will be confusion about +% which line is which. To fix this problem, one would want to be able to +% pick a much larger set of distinct colors, where the number of colors +% equals or exceeds the number of lines you want to plot. Because our +% ability to distinguish among colors has limits, one should choose these +% colors to be "maximally perceptually distinguishable." +% +% This function generates a set of colors which are distinguishable +% by reference to the "Lab" color space, which more closely matches +% human color perception than RGB. Given an initial large list of possible +% colors, it iteratively chooses the entry in the list that is farthest (in +% Lab space) from all previously-chosen entries. While this "greedy" +% algorithm does not yield a global maximum, it is simple and efficient. +% Moreover, the sequence of colors is consistent no matter how many you +% request, which facilitates the users' ability to learn the color order +% and avoids major changes in the appearance of plots when adding or +% removing lines. +% +% Syntax: +% colors = distinguishable_colors(n_colors) +% Specify the number of colors you want as a scalar, n_colors. This will +% generate an n_colors-by-3 matrix, each row representing an RGB +% color triple. If you don't precisely know how many you will need in +% advance, there is no harm (other than execution time) in specifying +% slightly more than you think you will need. +% +% colors = distinguishable_colors(n_colors,bg) +% This syntax allows you to specify the background color, to make sure that +% your colors are also distinguishable from the background. Default value +% is white. bg may be specified as an RGB triple or as one of the standard +% "ColorSpec" strings. You can even specify multiple colors: +% bg = {'w','k'} +% or +% bg = [1 1 1; 0 0 0] +% will only produce colors that are distinguishable from both white and +% black. +% +% colors = distinguishable_colors(n_colors,bg,rgb2labfunc) +% By default, distinguishable_colors uses the image processing toolbox's +% color conversion functions makecform and applycform. Alternatively, you +% can supply your own color conversion function. +% +% Example: +% c = distinguishable_colors(25); +% figure +% image(reshape(c,[1 size(c)])) +% +% Example using the file exchange's 'colorspace': +% func = @(x) colorspace('RGB->Lab',x); +% c = distinguishable_colors(25,'w',func); + +% Copyright 2010-2011 by Timothy E. Holy + + % Parse the inputs + if (nargin < 2) + bg = [1 1 1]; % default white background + else + if iscell(bg) + % User specified a list of colors as a cell aray + bgc = bg; + for i = 1:length(bgc) + bgc{i} = parsecolor(bgc{i}); + end + bg = cat(1,bgc{:}); + else + % User specified a numeric array of colors (n-by-3) + bg = parsecolor(bg); + end + end + + % Generate a sizable number of RGB triples. This represents our space of + % possible choices. By starting in RGB space, we ensure that all of the + % colors can be generated by the monitor. + n_grid = 30; % number of grid divisions along each axis in RGB space + x = linspace(0,1,n_grid); + [R,G,B] = ndgrid(x,x,x); + rgb = [R(:) G(:) B(:)]; + if (n_colors > size(rgb,1)/3) + error('You can''t readily distinguish that many colors'); + end + + % Convert to Lab color space, which more closely represents human + % perception + if (nargin > 2) + lab = func(rgb); + bglab = func(bg); + else + C = makecform('srgb2lab'); + lab = applycform(rgb,C); + bglab = applycform(bg,C); + end + + % If the user specified multiple background colors, compute distances + % from the candidate colors to the background colors + mindist2 = inf(size(rgb,1),1); + for i = 1:size(bglab,1)-1 + dX = bsxfun(@minus,lab,bglab(i,:)); % displacement all colors from bg + dist2 = sum(dX.^2,2); % square distance + mindist2 = min(dist2,mindist2); % dist2 to closest previously-chosen color + end + + % Iteratively pick the color that maximizes the distance to the nearest + % already-picked color + colors = zeros(n_colors,3); + lastlab = bglab(end,:); % initialize by making the "previous" color equal to background + for i = 1:n_colors + dX = bsxfun(@minus,lab,lastlab); % displacement of last from all colors on list + dist2 = sum(dX.^2,2); % square distance + mindist2 = min(dist2,mindist2); % dist2 to closest previously-chosen color + [~,index] = max(mindist2); % find the entry farthest from all previously-chosen colors + colors(i,:) = rgb(index,:); % save for output + lastlab = lab(index,:); % prepare for next iteration + end +end + +function c = parsecolor(s) + if ischar(s) + c = colorstr2rgb(s); + elseif isnumeric(s) && size(s,2) == 3 + c = s; + else + error('MATLAB:InvalidColorSpec','Color specification cannot be parsed.'); + end +end + +function c = colorstr2rgb(c) + % Convert a color string to an RGB value. + % This is cribbed from Matlab's whitebg function. + % Why don't they make this a stand-alone function? + rgbspec = [1 0 0;0 1 0;0 0 1;1 1 1;0 1 1;1 0 1;1 1 0;0 0 0]; + cspec = 'rgbwcmyk'; + k = find(cspec==c(1)); + if isempty(k) + error('MATLAB:InvalidColorString','Unknown color string.'); + end + if k~=3 || length(c)==1, + c = rgbspec(k,:); + elseif length(c)>2, + if strcmpi(c(1:3),'bla') + c = [0 0 0]; + elseif strcmpi(c(1:3),'blu') + c = [0 0 1]; + else + error('MATLAB:UnknownColorString', 'Unknown color string.'); + end + end +end diff --git a/eig2.m b/eig2.m new file mode 100644 index 0000000..90cf1f5 --- /dev/null +++ b/eig2.m @@ -0,0 +1,34 @@ +function D = eig2(A) +% function D = eig2(A) +% +% Compute in one shot the eigen-values of multiples (2 x 2) matrices +% +% INPUT: +% A: (2 x 2 x n) array +% OUTPUT: +% D: (2 x n). EIG2 returns in D(:,k) three eigen-values of A(:,:,k) +% +% See also: ParabolaRoots, eig3, eig +% +% Author: Bruno Luong +% History: +% Original 27-May-2010 + +if size(A,1) ~= 2 || size(A,2) ~= 2 + error('A must be [3x3xn] array'); +end + +A = reshape(A, 4, []).'; + +P3 = 1; +% Trace +P2 = -(A(:,1)+A(:,4)); + +% Determinant +P1 = A(:,1).*A(:,4) - A(:,2).*A(:,3); + +% Find the roots of characteristic polynomials +D = ParabolaRoots(P3, P2, P1).'; + +end % eig2 + diff --git a/ellipse.m b/ellipse.m new file mode 100644 index 0000000..202e6cd --- /dev/null +++ b/ellipse.m @@ -0,0 +1,149 @@ +function h=ellipse(ra,rb,ang,x0,y0,C,Nb) +% Ellipse adds ellipses to the current plot +% +% ELLIPSE(ra,rb,ang,x0,y0) adds an ellipse with semimajor axis of ra, +% a semimajor axis of radius rb, a semimajor axis of ang, centered at +% the point x0,y0. +% +% The length of ra, rb, and ang should be the same. +% If ra is a vector of length L and x0,y0 scalars, L ellipses +% are added at point x0,y0. +% If ra is a scalar and x0,y0 vectors of length M, M ellipse are with the same +% radii are added at the points x0,y0. +% If ra, x0, y0 are vectors of the same length L=M, M ellipses are added. +% If ra is a vector of length L and x0, y0 are vectors of length +% M~=L, L*M ellipses are added, at each point x0,y0, L ellipses of radius ra. +% +% ELLIPSE(ra,rb,ang,x0,y0,C) +% adds ellipses of color C. C may be a string ('r','b',...) or the RGB value. +% If no color is specified, it makes automatic use of the colors specified by +% the axes ColorOrder property. For several circles C may be a vector. +% +% ELLIPSE(ra,rb,ang,x0,y0,C,Nb), Nb specifies the number of points +% used to draw the ellipse. The default value is 300. Nb may be used +% for each ellipse individually. +% +% h=ELLIPSE(...) returns the handles to the ellipses. +% +% as a sample of how ellipse works, the following produces a red ellipse +% tipped up at a 45 deg axis from the x axis +% ellipse(1,2,pi/8,1,1,'r') +% +% note that if ra=rb, ELLIPSE plots a circle +% + +% written by D.G. Long, Brigham Young University, based on the +% CIRCLES.m original +% written by Peter Blattner, Institute of Microtechnology, University of +% Neuchatel, Switzerland, blattner@imt.unine.ch + + +% Check the number of input arguments + +if nargin<1, + ra=[]; +end; +if nargin<2, + rb=[]; +end; +if nargin<3, + ang=[]; +end; + +%if nargin==1, +% error('Not enough arguments'); +%end; + +if nargin<5, + x0=[]; + y0=[]; +end; + +if nargin<6, + C=[]; +end + +if nargin<7, + Nb=[]; +end + +% set up the default values + +if isempty(ra),ra=1;end; +if isempty(rb),rb=1;end; +if isempty(ang),ang=0;end; +if isempty(x0),x0=0;end; +if isempty(y0),y0=0;end; +if isempty(Nb),Nb=300;end; +if isempty(C),C=get(gca,'colororder');end; + +% work on the variable sizes + +x0=x0(:); +y0=y0(:); +ra=ra(:); +rb=rb(:); +ang=ang(:); +Nb=Nb(:); + +if isstr(C),C=C(:);end; + +if length(ra)~=length(rb), + error('length(ra)~=length(rb)'); +end; +if length(x0)~=length(y0), + error('length(x0)~=length(y0)'); +end; + +% how many inscribed elllipses are plotted + +if length(ra)~=length(x0) + maxk=length(ra)*length(x0); +else + maxk=length(ra); +end; + +% drawing loop + +for k=1:maxk + + if length(x0)==1 + xpos=x0; + ypos=y0; + radm=ra(k); + radn=rb(k); + if length(ang)==1 + an=ang; + else + an=ang(k); + end; + elseif length(ra)==1 + xpos=x0(k); + ypos=y0(k); + radm=ra; + radn=rb; + an=ang; + elseif length(x0)==length(ra) + xpos=x0(k); + ypos=y0(k); + radm=ra(k); + radn=rb(k); + an=ang(k); + else + rada=ra(fix((k-1)/size(x0,1))+1); + radb=rb(fix((k-1)/size(x0,1))+1); + an=ang(fix((k-1)/size(x0,1))+1); + xpos=x0(rem(k-1,size(x0,1))+1); + ypos=y0(rem(k-1,size(y0,1))+1); + end; + + co=cos(an); + si=sin(an); + the=linspace(0,2*pi,Nb(rem(k-1,size(Nb,1))+1,:)+1); +% x=radm*cos(the)*co-si*radn*sin(the)+xpos; +% y=radm*cos(the)*si+co*radn*sin(the)+ypos; + h(k)=line(radm*cos(the)*co-si*radn*sin(the)+xpos,radm*cos(the)*si+co*radn*sin(the)+ypos); + set(h(k),'color',C(rem(k-1,size(C,1))+1,:),'linewidth',2,'linestyle',':'); + +end; + diff --git a/finding_true_num_2D.m b/finding_true_num_2D.m new file mode 100644 index 0000000..27afdc1 --- /dev/null +++ b/finding_true_num_2D.m @@ -0,0 +1,19 @@ +function [dis_min, true_num] = finding_true_num_2D(p0,p_s) + +% ========================================================================= +% sub function called by iterative particle reconstruction +% determine the true particle in the original volume +% that the reconstructed particle corresponds to. +% by Tianqi Guo, Aug 2017 +% ========================================================================= + +% calculate the distance in x, y, z the reconstructed particle with +% all the true particles in the original volume +du_i = (p0(:,1) - p_s(ones(size(p0,1),1),1)); +dv_i = (p0(:,2) - p_s(ones(size(p0,1),1),2)); + +% finding the closest true particle and return the distance +g_t = sum([du_i dv_i].^2,2).^0.5; +[dis_min, true_num] = min(g_t); + +end \ No newline at end of file diff --git a/parsave.m b/parsave.m new file mode 100644 index 0000000..b3c18d6 --- /dev/null +++ b/parsave.m @@ -0,0 +1,3 @@ +function parsave(path,data) +save(path,'data') +end \ No newline at end of file diff --git a/polartrans.m b/polartrans.m new file mode 100644 index 0000000..cfc621c --- /dev/null +++ b/polartrans.m @@ -0,0 +1,89 @@ +% POLARTRANS - Transforms image to polar coordinates +% +% Usage: pim = polartrans(im, nrad, ntheta, cx, cy, linlog, shape) +% +% Arguments: +% im - image to be transformed. +% nrad - number of radius values. +% ntheta - number of theta values. +% cx, cy - optional specification of origin. If this is not +% specified it defaults to the centre of the image. +% linlog - optional string 'linear' or 'log' to obtain a +% transformation with linear or logarithmic radius +% values. linear is the default. +% shape - optional string 'full' or 'valid' +% 'full' results in the full polar transform being +% returned (the circle that fully encloses the original +% image). This is the default. +% 'valid' returns the polar transform of the largest +% circle that can fit within the image. +% +% Returns pim - image in polar coordinates with radius increasing +% down the rows and theta along the columns. The size +% of the image is nrad x ntheta. Note that theta is +% +ve clockwise as x is considered +ve along the +% columns and y +ve down the rows. +% +% When specifying the origin it is assumed that the top left pixel has +% coordinates (1,1). + +% Copyright (c) 2002 Peter Kovesi +% School of Computer Science & Software Engineering +% The University of Western Australia +% http://www.csse.uwa.edu.au/ +% +% Permission is hereby granted, free of charge, to any person obtaining a copy +% of this software and associated documentation files (the "Software"), to deal +% in the Software without restriction, subject to the following conditions: +% +% The above copyright notice and this permission notice shall be included in +% all copies or substantial portions of the Software. +% +% The Software is provided "as is", without warranty of any kind. + +% December 2002 +% November 2006 Correction to calculation of maxlogr (thanks to Chang Lei) + +function [theta, radius, pim] = polartrans(im, nrad, ntheta, cx, cy, linlog, shape) + +[rows, cols] = size(im); + +if nargin==3 % Set origin to centre. + cx = cols/2+.5; % Add 0.5 because indexing starts at 1 + cy = rows/2+.5; +end + +if nargin < 7, shape = 'full'; end +if nargin < 6, linlog = 'linear'; end + +if strcmp(shape,'full') % Find maximum radius value + dx = max([cx-1, cols-cx]); + dy = max([cy-1, rows-cy]); + rmax = sqrt(dx^2+dy^2); +elseif strcmp(shape,'valid') % Find minimum radius value + rmax = min([cx-1, cols-cx, cy-1, rows-cy]); +else + error('Invalid shape specification'); +end + +% Increments in radius and theta + +deltatheta = 2*pi/ntheta; + +if strcmp(linlog,'linear') + deltarad = rmax/(nrad-1); + [theta, radius] = meshgrid([0:ntheta-1]*deltatheta, [0:nrad-1]*deltarad); +elseif strcmp(linlog,'log') + maxlogr = log(rmax); + deltalogr = maxlogr/(nrad-1); + [theta, radius] = meshgrid([0:ntheta-1]*deltatheta, exp([0:nrad-1]*deltalogr)); +else + error('Invalid radial transformtion (must be linear or log)'); +end + +xi = radius.*cos(theta) + cx; % Locations in image to interpolate data +yi = radius.*sin(theta) + cy; % from. + +[x,y] = meshgrid([1:cols],[1:rows]); +pim = interp2(x, y, double(im), xi, yi); + diff --git a/socdiff.m b/socdiff.m new file mode 100644 index 0000000..f05c359 --- /dev/null +++ b/socdiff.m @@ -0,0 +1,56 @@ +function [N]=socdiff(N,dx,dir) + +size_N = size(N); + +if isempty(dir) + if size_N(1)~=1 + dir = 1; + elseif size_N(2)~=1 + dir = 2; + elseif size_N(3)~=1 + dir = 3; + else + dir = 1; + end +end + +ndim_N = ndims(N); +if ndim_N > 3 + error('function only defined up to 3D matrices') +end + +if isempty(dx) + dx = 1; +end + + +N = permute(N, [dir, 1:dir-1, dir+1:ndim_N]); +size_N = size(N); %find again for permuted matrix +ndim_N = ndims(N); + + +NI = size_N(1); +NJ = size_N(2); +if ndim_N == 3 + NK = size_N(3); +else + NK = 1; +end + + +if NI < 3 + error('stencil size is 3 pts in active direction') +end + +DN = zeros(NI,NJ,NK); +%for j=1:NJ +% for k=1:NK + DN(1 ,:,:) = ( -N( 3 ,:,:) +4*N( 2 ,:,:) -3*N(1 ,:,:)) / (2*dx); + DN(2:NI-1,:,:) = ( N( 3:NI,:,:) - N(1:NI-2,:,:)) / (2*dx); + DN(NI ,:,:) = (3*N(NI ,:,:) -4*N(NI-1,:,:) + N( NI-2,:,:)) / (2*dx); +% end +%end + +%keyboard + +N = ipermute(DN, [dir, 1:dir-1, dir+1:ndim_N]); diff --git a/weighted_nearest_neighbor3D_MP.m b/weighted_nearest_neighbor3D_MP.m new file mode 100644 index 0000000..03fa20e --- /dev/null +++ b/weighted_nearest_neighbor3D_MP.m @@ -0,0 +1,137 @@ +function tracks = weighted_nearest_neighbor3D_MP(p2_est,p1_pos,p2_pos,paras_1,paras_2,s_radius,weights) + +% modified based on the multi-parametric PTV by Cardwell's paper +% tracking objects defined by +% 1) position (x,y,z) that changes over frames +% 2) property parameters (p_1,2,3,...) that remain almost unchanged between +% adjacent frames. +% This function finds the most probable matchings between the tracking +% objects from two frames by looking at the minimum of weighted +% difference calculated from position and property info. + +% inputs: +% 1) p2_est, estimated positions of objects in the 2nd frame +% 2) p1_pos,p2_pos, positions of objects in the two frames +% 3) paras_1,paras_2, properties of objects in the two frames +% 4) s_radius, search radius from the estimated position +% 5) weights, relative importance of each property +% outputs: +% 1) tracks, the most probable matching of objects between two frames. + +% this code is part of Prana from Vlachos reseach group at Purdue +% University +% by Tianqi Guo. (Jan. 2018) +%========================================================================== + +% number of tracking objects from two frames +num_p1 = length(p1_pos); num_p2 = length(p2_pos); + +% number of property parameters +num_paras = size(paras_1,2); + +% maximum difference of each property parameter +paras_diff_max = zeros(1,num_paras); + +% calculate the maximum difference for each parameter +for i = 1 : num_paras + paras_diff_max(i) = max([paras_1(:,i);paras_2(:,i)]) - min([paras_1(:,i);paras_2(:,i)]); +end + +% if the parameter is the same for all objects, +% exclude it from the difference calculation +paras_diff_max(paras_diff_max == 0) = Inf; + +% record the actual positions +p1_org = p1_pos; +p1_pos = p2_est; + +% placeholder +compare = cell(num_p1,1); + +% for each object in the 1st frame +% could be replaced with regular for loop +for i = 1:num_p1 + + % difference in positions\distance between this object and all the + % objects in the 2nd frame + dX = p2_pos(:,1) - p1_pos(i,1); + dY = p2_pos(:,2) - p1_pos(i,2); + dZ = p2_pos(:,3) - p1_pos(i,3); + distance = sqrt(dX.^2+dY.^2+dZ.^2); + + % only consider the objects within the search radius + p_index = find(distance<=s_radius); + compare_i = zeros(length(p_index),1); + % structure of compare_i + % [object # in 1st frame,object # in 2nd frame,weighted difference] + + % This statement catches the possibility that no paritlces are with in + % the search radius. + if ~isempty(p_index) + + % object # in 1st frame + compare_i(:,1) = ones(size(p_index)).*i; + + % object # in 2nd frame + compare_i(:,2) = p_index; + + % matching probablility\weighted difference + Prob = zeros(size(p_index,1),1+num_paras); + + % compute the match probability\weighted difference for each possible pairing + Prob(:,1) = (distance(p_index)./s_radius).*weights(1); + % for each tracking parameter + for j = 1:num_paras + + % normalized by the maximum difference of this parameter + Prob(:,j+1) = (abs(paras_2(p_index,j) - paras_1(i,j))./paras_diff_max(j)).*weights(j+1); + + end + + % the weighted difference + compare_i(:,3)=sum(Prob,2)./sum(weights); + + % populate main 'compare' array + compare{i}=compare_i; + end + +% clear dX dY dZ distance p_index compare_i Prob +end +clear i +compare = cell2mat(compare); + +% Check to make sure that the compare variable is not empty. It will be +% empty if no particles where paired durring this processes. The use +% should try increasing there search radius. +if ~isempty(compare) + %sort the comparison array in ascending order based on the pairing + %probability + compare_sort=sortrows(compare,3); + + %determine the best particle matches and successively match particles + %until the max number of pairs has been reached (or no more exist) + p_pairs=[]; + max_matches=min(num_p1,num_p2); + found_matches_im1=[]; found_matches_im2=[]; + c=0; i=1; + while (c<=max_matches) && (i<=size(compare_sort,1)) + test1=found_matches_im1==compare_sort(i,1); + test2=found_matches_im2==compare_sort(i,2); + if (any(test1) || any(test2)) + i=i+1; + else + found_matches_im1=vertcat(found_matches_im1,compare_sort(i,1)); + found_matches_im2=vertcat(found_matches_im2,compare_sort(i,2)); + p_pairs=vertcat(p_pairs,compare_sort(i,:)); + i=i+1; + c=c+1; + end + end + + %populate the 'tracks' arrays with the following structure: + % tracks=[X1, X2, Y1, Y2, Z1, Z2, p#1 p#2 match_probability] + tracks = [p_pairs(:,:), p1_org(p_pairs(:,1),:), p2_pos(p_pairs(:,2),:)]; +else + tracks = []; +end +