format compact clearvars, clc close all warning('off','all') set(0,'DefaultFigureWindowStyle','normal') set(0,'DefaultAxesTitleFontWeight','normal'); set(0,'DefaultAxesFontName','Avenir') %% this data table defines the key settings Sguv = readtable('guv_settings.csv'); runtype = 2; % 1: run just testing few frames % 2: run full kymograph analysis on 1 or more GUVs % 3: find and load previous data and just adjust plots % 4: manually select file and adjust % name of output folder - ALWAYS HAVE BACKSLASH datatype = 'NOlumen_60plperp_SAVEALL_v2022-11-22/'; % choose which guvs to analyze Gidx_analysis = [1:18];%1,5,9,12,13,15:17 n_analysis = length(Gidx_analysis); % # of GUVs to be analyzed savefigs = 2; % save figures and movies % 0 = no; % 1= mfiles,kymos,figs ; % 2 = 1 + dataRaw(images,other boundary quants etc) %% if runtype == 1 % testing, no kymos % guvname = 'G01'; % listFrames = 1:28; % t_rapa = 7; % nsize = 1; % thetacenter = 0; % only single guv in gidx_analysis for testing ii = Gidx_analysis(1); guvname = Sguv.guvname{ii}; lastframe = Sguv.lastframe(ii); listFrames = 1:lastframe; t_rapa = Sguv.t_rapa(ii); nsize = Sguv.nsize(ii); thetacenter = Sguv.thetacenter(ii); runall(runtype,datatype,guvname,t_rapa,thetacenter,listFrames,nsize,savefigs); elseif runtype == 2 % MAIN RUNS for i = 1:n_analysis close all; ii = Gidx_analysis(i); guvname = Sguv.guvname{ii}; lastframe = Sguv.lastframe(ii); listFrames = 1:lastframe; t_rapa = Sguv.t_rapa(ii); nsize = Sguv.nsize(ii); thetacenter = Sguv.thetacenter(ii); runall(runtype,datatype,guvname,t_rapa,thetacenter,listFrames,nsize,savefigs); end elseif runtype == 3 % load previously saved data % list of ALL folder Sfolds = dir('G*output'); % list of Global output folders df_list = cell(n_analysis,1); % list of .mat files for SELECT GUVs %df = cell(n_analysis,1); % variables in each .mat file tstampnew = datestr(now,'dd_mmm_yyyy_HH_MM_SS'); for i = 1:n_analysis close all; tmpfolder = fullfile(Sfolds(Gidx_analysis(i)).name,datatype); tmpfile = dir([tmpfolder,'*.mat']); df_list{i} = fullfile(tmpfolder,tmpfile.name); load(df_list{i},'dataRaw','dataParams','dataKymoBin','dataKymoVar'); % adjust this to ensure internal saving of png within function goes to correct folder dataParams.outputfolder = tmpfolder; kf = dataParams.kymofixed; %% add mask timeseries % input variables tstampnewtext = strcat(tstampnew,'_MASKTS'); dataParams.tstamp1 = tstampnewtext; tmpKymo = dataKymoBin; makekymo = dataParams.makekymo; % plot movies and save only based on the first kymograph set (fixedpts or fixedtheta usually) plottseries(dataParams,dataRaw,makekymo,tmpKymo,kf(1)); %% ADJUST FONT SIZE % tstampnewtext = strcat(tstampnew,'_FONT14'); % dataParams.tstamp1 = tstampnewtext; % % count = 0; % png saving step depends on count % adjfix(1) = plotphyskymos(dataParams,dataKymoBin,kf(1),count); % adjfix(2) = plotbiochemkymos(dataParams,dataKymoBin,kf(1),count); % count = 5; % adjvar(1) = plotphyskymos(dataParams,dataKymoVar,kf(2),count); % adjvar(2) = plotbiochemkymos(dataParams,dataKymoVar,kf(2),count); % if savefigs > 0 % %savefig(f01,strcat(dataParams.outputfolder,dataParams.guvname,'_',tstampnewtext,'_NEW_macrosummary.fig'),'compact'); % savefig(adjfix,strcat(dataParams.outputfolder,dataParams.guvname,'_',tstampnewtext,'_figsoutFix.fig'),'compact'); % savefig(adjvar,strcat(dataParams.outputfolder,dataParams.guvname,'_',tstampnewtext,'_figsoutVar.fig'),'compact'); % end %% ADJUST BINARIZATION FOR PATTERN VISUALIZATION % % adjust time stamp to include adjustment % tstampnewtext = strcat(tstampnew,'_PATTERNS'); % dataParams.tstamp1 = tstampnewtext; % % count = 10; % % kmemb = dataKymoBin.kymofluor{1}; % kacta = dataKymoBin.kymofluor{3}; % acta % kactin = dataKymoBin.kymofluor{2}; % actin % kcurv = dataKymoBin.kymocurv; % % actinthresh = median(kactin(:)); % % scrv = sort(abs(kcurv(:)),'descend','MissingPlacement','last'); % ncrv = length(scrv); % crvthresh = scrv(round((0.10)*ncrv)); % top 20% % % kactin1 = (kactin>=actinthresh); % threshold value different for each GUV % kcurv1 = (abs(kcurv)>=crvthresh); % %figure; imagesc(kactin1); % % %dataKymoBin.kymofluor{2} = kactin1; % dataKymoBin.kymocurv = kcurv1; % % % figsoutFix(1) = plotphyskymos(dataParams,dataKymoBin,kf(1),count); % % clrs = [0 0.4470 0.7410; 0.9290 0.6940 0.1250]; % % ax = gca; colorbar off % % ax.CLim = [0,1]; % % colormap(ax,clrs); % % colorbar % % % figsoutFix(2) = plotbiochemkymos(dataParams,dataKymoBin,kf(1),count); % % ax = gca; colorbar off % % ax.CLim = [0,1]; % % colormap(ax,clrs); % % colorbar % % count = 20; % % kactinvar = dataKymoVar.kymofluor{2}; % actin % kcurvvar = dataKymoVar.kymocurv; % % kactinvar1 = (kactinvar>=actinthresh)&(~isnan(kactinvar)); % threshold value different for each GUV % kcurvvar1 = (abs(kcurvvar)>=crvthresh)&(~isnan(kcurvvar)); % %figure; imagesc(kactin1); % % kactinVnull = isnan(kactinvar); % kactinVfin = double(kactinvar1); % kactinVfin(kactinVnull) = nan; % % kcurvVnull = isnan(kcurvvar); % kcurvVfin = double(kcurvvar1); % kcurvVfin(kcurvVnull) = nan; % % % % dataKymoVar.kymofluor{2} = kactinVfin; % dataKymoVar.kymocurv = kcurvVfin; % % figsoutVar(1) = plotphyskymos(dataParams,dataKymoVar,kf(2),count); % % ax = gca; colorbar off % % ax.CLim = [0,1]; % % colormap(ax,clrs); % % colorbar % % figsoutVar(2) = plotbiochemkymos(dataParams,dataKymoVar,kf(2),count); % % ax = gca; colorbar off % % ax.CLim = [0,1]; % % colormap(ax,clrs); % % colorbar % if savefigs > 0 % %savefig(f01,strcat(dataParams.outputfolder,dataParams.guvname,'_',tstampnewtext,'_NEW_macrosummary.fig'),'compact'); % savefig(figsoutFix,strcat(dataParams.outputfolder,dataParams.guvname,'_',tstampnewtext,'_figsoutFix.fig'),'compact'); % savefig(figsoutVar,strcat(dataParams.outputfolder,dataParams.guvname,'_',tstampnewtext,'_figsoutVar.fig'),'compact'); % end end elseif runtype == 4 % chose manually % close all; % tstampnew = datestr(now,'dd_mmm_yyyy_HH_MM_SS'); % % [filelist,filepath]=uigetfile('*.mat','MultiSelect','on',... % 'Select LIST to plot'); pause(0.5); % % tmpfullfile = fullfile(filepath,filelist); % load(tmpfullfile,'dataRaw','dataParams'); % % dataParams.do_fixedbox = 1; % % % adjust this to ensure internal saving of png within function goes to correct folder % dataParams.outputfolder = filepath; % dataParams.tstamp1 = strcat(tstampnew,'_NEW'); % % % f01 = plotmovement(dataParams,dataRaw); % % if savefigs == 2 % savefig(f01,strcat(dataParams.outputfolder,dataParams.guvname,'_',tstampnew,'_NEW_macrosummary.fig'),'compact'); % elseif savefigs == 1 % savefig(f01,strcat(dataParams.outputfolder,dataParams.guvname,'_',tstampnew,'_NEW_macrosummary.fig'),'compact'); % end end function runall(runtype,datatype,guvname,t_rapa,thetacenter,listFrames,nsize,savefigs) %-------------------------------------------------------------------------- % SET ITYPE IN INTERPBOUNDARY! 'pchip' or 'makima' % SET HOW FLUOR INTENSITIES ARE COMPUTED! % SET kymograph smoothing window % guvname = 'G09'; % t_rapa = 7; % last frame before rapamycin added, approx = time 0 for recruitment % thetacenter = 0; % angle (degrees rel. to x-axis) reference point to align center of kymograph % listFrames = [1:74]; % choose frames (not necessarily consecutive) to include in analysis, in ascending order % %IMPORTANT FOR SETTING SIZE OF INNER OUTER MASKS!!! % if strcmp(guvname(1),"G"); nsize = 0.7; % 1 for G01, % else; nsize = 3; % 3 for L48, check tseries inner/outer masks % end % m0 = 3; % boundary points per pixel; will multiply perimeter, used in pre-computation interp m1 = 6; % muliple of perimeter, post-computation interp - needed for small guvs to get 301 spatial bins % basystem = 'mac'; kymofixed = [2 0]; % 0: varying points; % 1: fixedpts: use fixed # of points for 2nd kymograph (= # pts in longest perim, post-interp) % 2: fixedtheta: % use fixed # of angles for 2nd kymograph nbins = 360; % ODD # of bins to use if kymofixedpts or kymofixedtheta pvary = 3; % multiple of perimeter for final interpolation of varying kymo % ODD ensures needle angle will be captured in 1 bin at center % 0/1 or true/false: do_plotAlignseries = 0; % display boundary, centroid, kymo center point do_plotTimeseries = 0; % display MAIN time series : FLUORESC do_plotMaskseries = 0; % plot inner, outer mask time series do_plotBoundaryOverlay = 0; % separate plot of boundary overlay over time % more defaults 0/1 here plperp = 0.6; % 60% outward from inner line for intensity computation do_top3px = 0; % use only top 3 px for fluor intensity do_membnorm = 0; % divide biochem signals by membrane channel do_fixedbox = 1; % in summary fig, plot with fixed box size for centroid and boundary evolution do_plotVelseries = 0; % very slow - plotting displacements with velocity in color do_plotEllipseries = 0; % plot bounding ellipse do_aspratio_ind = 1; % 1 = do independent asp ratio computation, % 0 = use matlab region props maj/min axis to compute aspect ratio veltrack1 = 1; % compute velocity using motion tracking (1) or pdist2 (0) for fixed kymo veltrack2 = 0; % compute velocity using motion tracking (1) or pdist2 (0) for variable kymo % alignvel will be based on whatever is chosen here. bdyvel is still pdist2 based % Select curvature calculation method computecurv = 2; % 1 for circumcenter; 2 for fitcircle % Select smoothing savgolay_smooth = true; % do savitzky-golay smoothing of boundary for curvature estimates darc = 15; % distance (pixels) of arc segment used for savgolay darc2 = 15; % distance (pixels) of arc segment used for curvature computation % parc = 0.05; % prop of perimeter: arc segment used for savgolay % parc2 = 0.05; % prop of perimeter: arc segment used for curvature computation do_movavg = 0; % do moving average s_movavg = 2*m1; % if do_movavg, use this for # of points right/left %------------------------------------------ minfont = 18; sat = 1; % factor for caxis upper limit fluor signals satlow = 0; % factor for caxis lower limit of fluor signal %-------------------------------- inputfolder = strcat(guvname,'_input/',guvname,'_acs/'); %--------------------- if runtype == 1 makekymo = 1; % create kymograph outputfolder = strcat(guvname,'_testout/',datatype); % name of output folder elseif runtype == 2 makekymo = 1; % Create kymograph outputfolder = strcat(guvname,'_output/',datatype); % name of output folder else % only re-runs makekymo = 0; % No kymograph outputfolder = strcat(guvname,'_output/',dataype); end if strcmp(guvname(1),'L') fname1 = strcat(guvname,'_mCh_memb'); % membrane marker file name fname2 = strcat(guvname,'_YFP_actin'); % actin marker file name fname3 = strcat(guvname,'_CFP_ActA'); fname4 = strcat(guvname,'_BF'); fname5 = strcat(guvname,'_647_dye'); fnames = {fname1;fname2;fname3;fname4;fname5}; needlech = 4 ; % BF channel for needle elseif strcmp(guvname(1),'G') fname1 = strcat(guvname,'_mCh_memb'); % membrane marker file name fname2 = strcat(guvname,'_YFP_actin'); % actin marker file name fname3 = strcat(guvname,'_CFP_ActA'); fnames = {fname1;fname2;fname3}; end maskfname = strcat(guvname,'_mask'); membch = 1 ; % channel of membrane marker fluorch = [1:3]; % membrane 1, actin 2, ActA 3; nfluor = length(fluorch); nchannels = length(fnames); fullFileNames = cell(nchannels,1); for i = 1:nchannels fullFileNames{i} = strcat(inputfolder,fnames{i},'_final.tif'); % if ~isfile(fullFileNames{i}) % fullFileNames{i} = strcat(inputfolder,fnames{i},'.tif'); % end end tstamp1 = datestr(now,'dd_mmm_yyyy_HH_MM_SS'); finfo = imfinfo(fullFileNames{1}); nallframes = length(finfo); cols = finfo.Width; rows = finfo.Height; bits = finfo(1).BitDepth; im0 = listFrames(1); imL = min(nallframes,listFrames(end)); deltaf = imL - im0 + 1; % useful when skipping frames during testing nframes = min(nallframes,length(listFrames)); dt = 1; % keep 1 for pixels/frame; use seconds per frame for pixels/sec % dx = 1/finfo.Xresolution; %1/2.2764 = 0.44 microns/pixel if exist(outputfolder,'dir') == 0 mkdir(outputfolder); end %--------------------------- %ticks for kymographs if nframes < 10 tickstep = 1; else tickstep = floor(nframes/10); end %-------------------------------------------------------------------------- maskFileName = strcat(inputfolder,maskfname,'_final.tif'); if ~isfile(maskFileName) maskFileName = strcat(inputfolder,maskfname,'.tif'); end maskInfo = imfinfo(maskFileName); if nallframes~=length(maskInfo) disp('Error: Number of frames in the image and mask files must be the same') end % save initial settings dataParams = v2struct; %---------------------------------------------------------- % cell arrays that are not channel specific allCurvatures = cell(nframes,1); % Curvatures, all frames allpdist2Velocities = cell(nframes,1); % Boundary velocity, all fr. allpdist2Dcumuls = cell(nframes,1); allBoundaries = cell(nframes,1); % Boundary positions, all fr. % channel specific allImages = cell(nframes,nchannels); % original images, frame i, channel j Igs = cell(nframes, nchannels); Iinit = cell(nframes, nchannels); allIrgb = cell(nframes, nchannels); allIntensities = cell(nframes,nfluor); % boundary intensities for all frames and fluorescent channels alliBW = cell(nframes,nfluor); alloBW = cell(nframes,nfluor); allinner = cell(nframes,nfluor); allouter = cell(nframes,nfluor); % numeric arrays allCentroids = zeros(nframes,2); % Centroid, all frames allAspRatio_ml = zeros(nframes,1); % compute using matlab ratio of major/minor axis length from region props allAspRatio_ind = zeros(nframes,1); % independent computation of Aspect ratio (ratio of major/minor axis of minimum bounding ellipse) allA0 = cell(nframes,1); % independent computation allows for plotting the best fit ellipse allA0c = cell(nframes,1); % independent computation allows for plotting the best fit ellipse allEccents = zeros(nframes,1); %Eccentricity allMajAxL = zeros(nframes,1); % Major axis length allLengths = zeros(nframes,1); % Length of boundaries allPerims = zeros(nframes,1); allAreas = zeros(nframes,1); % Cell areas, all frames allBW = zeros(rows,cols,nframes,'logical'); % Segmentation all frames - only based on 1 channel %% compute boundary points and smooth %---------------------------------------------------- % This section does velocity/boundary computations for membrane channel % and intensity computations for every channel for i = 1:nframes movieFrame=listFrames(i); %im0:imL; % TIFF indices frame = i; %movieFrame-im0+1; % indices for saving data fprintf('Frame number: %4d. ',movieFrame); for k = 1:nchannels %--------------------------------------------- % gaussian smoothing for membrane (I) and actin (I2) [Xgs,Xorig] = readframe(fullFileNames{k},movieFrame); %uint16 output (not rgb) allImages{frame,k} = Xorig; Igs{frame,k} = Xgs; % define initial images for computations % segment based on I (gaussian smoothed membrane maker) % or segment based on Iorig (original membrane marker) Xinit = Xgs; Iinit{frame,k} = Xinit; % make boundary marker into rgb image for plotting % use imadjust to increase contrast Xrgb = repmat(imadjust(Xinit),1,1,3); allIrgb{frame,k} = Xrgb; end BW = logical(imread(maskFileName,movieFrame)); [B,L,ncells] = bwboundaries(BW,'noholes'); stats = regionprops(L,'Area','Centroid','Perimeter','Eccentricity','MajorAxisLength','MinorAxisLength'); % identify boundary points for single object (exclude new ones) if ncells==1 % BOUNDARY FROM MASK bdyraw = B{1}; % (row,col) coordinates (y,x) area = stats.Area; centroid = stats.Centroid; perimeter = stats.Perimeter; eccent = stats.Eccentricity; mjaxl = stats.MajorAxisLength; aspratio_ml = stats.MajorAxisLength/stats.MinorAxisLength; % interpolate boundary to have spacing of 1 pixel bdy = interpboundary(bdyraw,m0*perimeter+1); elseif ncells > 1 && frame > 1 tmpcenters = zeros(ncells,2); tmpdist = zeros(ncells,1); for k=1:ncells tmpcenters(k,:) = stats(k).Centroid; tmpdist(k) = norm(tmpcenters(k,:)-allCentroids(frame-1,:)); end [~,kstar] = min(tmpdist); % BOUNDARY FROM MASK bdyraw = B{kstar}; % (row,col) coordinates (y,x) area = stats(kstar).Area; centroid = stats(kstar).Centroid; perimeter = stats(kstar).Perimeter; eccent = stats(kstar).Eccentricity; mjaxl = stats(kstar).MajorAxisLength; aspratio_ml = stats(kstar).MajorAxisLength/stats(kstar).MinorAxisLength; % interpolate boundary to have spacing of 1 pixel; bdy = interpboundary(bdyraw,m0*perimeter+1); % #o unique points + 1 end %--------------------------------------------- % optional savitzky-golay smoothing of boundary % Now smooth with a Savitzky-Golay sliding polynomial filter %deltaarc1 = parc*perimeter*m0; % #points in arc = perimeter fraction * points per pixel (m0) deltaarc1 = darc*m0; % #points in arc = arc distance (pixels) * points per pixel (m0) windowWidth = 2*(round(0.5*deltaarc1))+1; % must be odd; adjust based on length of interpBdy polynomialOrder = 2; % keep this as 2 for quadratic smoothing of boundary % use membrane channel rgb to display bdy_sg = sgsmoothboundary(allIrgb{frame,membch},bdy,windowWidth,polynomialOrder); %-------------------------------------- if movieFrame > im0 oBorigs = Borigs; % save old boundary for velocity before updating end %---------------------------------------- % ============ SELECT boundary to use if savgolay_smooth Borigs = bdy_sg; % savitzky-golay smoothed boundary % sgolay smoothed boundary will only be approximately circular % (not closed loop usually due to edge effects) else Borigs = bdy; % unsmoothed boundary end %---------------------------------------- %% compute aspect ratio %prm0 = v2struct(do_plotEllipseries,movieFrame,minfont); if do_aspratio_ind [aspratio_ind,A0,A0c,~] = aspectratio(Borigs); end %% needle position check %% compute velocities (pixels/frame) if movieFrame > im0 %&& ~veltrack % Get displacement between old boundary and new boundary over 1 frame % (row,col) coordinates (y,x) % dt = 1 so it's same as displacements pVel = v2struct(do_plotVelseries,savefigs,movieFrame,im0,outputfolder,... guvname,tstamp1); [vel,idx] = getvelocities(Borigs,oBorigs,dt,pVel,allIrgb{frame,membch}); dcumul = dcumul(idx) + dt*vel; %velold = vel; else vel = nan(length(Borigs),1); %velold = zeros(length(Borigs),1); dcumul = zeros(length(Borigs),1); end %% compute intensity from all channels % At the ith frame Iinit_curr = Iinit(frame,:); % cell array of all channels at current frame pFluor = v2struct(centroid,movieFrame,im0,imL,do_plotTimeseries,do_membnorm,do_top3px,savefigs,outputfolder,... tstamp1,fnames,fluorch,nfluor,nchannels,nsize,plperp,bits); [fluor,iBW,oBW,inner,outer] = getIntensity(Borigs,Iinit_curr,BW,pFluor); % get actin intensity +/- boundary % --------------- %% compute k-circumcenter % % =========== CIRCUMCENTER % % circumcenter approach to estimate curvature % % % define spaced boundary points for curvature estimation % %defBk = 'distance'; % define distance (in pixels) separating boundary points % defBk = 'number'; % define number of boundary points % %defBk= 'original'; % use original boundary points % % Bk_dist = parc2*perimeter*m0; %#points spaced out in arc = perimeter fraction * points per pixel (m0) % Bk_num = round(perimeter/deltaarc1); % number of unique (non-circular) spaced points around boundary % curvatureThresh = 1; %0.06 the maximum allowed value of the curvature measure % loopclose = 1; % 0 - if open boundaries | 1 - if closed boundaries % iflag = 2; % two curvature interpolation options for circumcenter approach % % needed b/c approach uses spaced boundary points % % 1 use scatteredinterpolant to match coordinates of bdy ; % % 2 use interpboundary to match length of bdy % % % spaced boundary points for estimating k (curvature) % % Bkpoints is circular % Bkpoints = defBkpoints(I1rgb,Borigs,defBk,Bk_dist,Bk_num); % % Pcurv1 = v2struct(curvatureThresh,loopclose,perimeter); % [Bkcurv0] = curvature1_v5(Bkpoints,Pcurv1); % % Bkcurv0 is NOT circular; Borigs,Bkpoints are circular % % make Bkcurv circular % Bkcurv = [Bkcurv0;Bkcurv0(1)]; % % % interpolate curvature to match length or coordinates of bdy (iflag = 1 or 2) % % and plot overlay as fig 21 % curv1 = plotcurvature1_v6(I1rgb,Borigs,Bkcurv,Bkpoints,iflag); % % ------------- %% compute k-fitcircle % ====================== FITCIRCLE % fitcircle approach to estimate curvature %deltaarc2 = parc2*perimeter*m0; % 0.1, %#points in arc = perimeter fraction * points per pixel (m0) deltaarc2 = darc2*m0; % 0.1, %#points in arc = arc distance (pixels) * points per pixel (m0) deltaarc2_half=round(0.5*deltaarc2); % 2*delta+1 is used for symmetry c2 = curve_fitcircle(Borigs,deltaarc2_half); % use savgolay smoothed or raw boundary % Borigs is the selected boundary for circumcenter/fitcircle = bdy or bdy_sg curv2 = plotcurvature2_v1(allIrgb{frame,membch},Borigs,c2); % plot on original boundary % test %c2 = curve_fitcircle(Bkpoints,1); % use savgolay smoothed or raw boundary %curv2 = plotcurvature2_v1(I1rgb,bdy,c2,Bkpoints,iflag); % Select which curvature approach to use if computecurv == 1 curv = curv1; % use circumcenter approach % show boundary spacing, smoothing used for circumcenter % f1 = displayboundary(I1rgb,bdy,Borigs,Bkpoints); else curv = curv2; % use fitcircle approach end %% post computation interpolation of all boundary measurements n = round(m1*perimeter)+1;% set n to be: (number of new UNIQUE boundary points) + 1 % if want spacing of 0.5 pixel unit, set # of unique boundary points = 2*perimeter % n = 2*perimeter + 1 % gfp, vel only approximately circular; bdy, curv are circ iFluor = cell(1,nfluor); % pre-allocation is required % interpolate bdy quantities % pass fluor as comma separated list, not as a cell array [iBdy,iVel,iCurv,iDcumul,iFluor{:}] = interpboundary(Borigs,n,vel,curv,dcumul,fluor{:}); % interpBdy and interpCurv are circular; others just very close % make circular all iVel(end) = iVel(1); iCurv(end) = iCurv(1); iDcumul(end)= iDcumul(1); matiFluor = cell2mat(iFluor); matiFluor(end,:) = matiFluor(1,:); iFluor = mat2cell(matiFluor,length(iFluor{1}),ones(1,nfluor)); if savefigs > 0 % save m file if movieFrame == imL % save mfile if saving figs currentfile = strcat(mfilename,'.m'); destinationfile = strcat(outputfolder,mfilename,'_',tstamp1,'.m'); copyfile(currentfile,destinationfile); end end if frame>1 allpdist2Velocities{frame} = iVel; % BSA else allpdist2Velocities{frame} = nan(length(iBdy),1); % BSA end allPerims(frame) = perimeter; allCurvatures{frame} = iCurv; allpdist2Dcumuls{frame} = iDcumul; allBoundaries{frame} = iBdy; allIntensities(frame,:) = iFluor; % row = frame, cols = nfluor allCentroids(frame,:) = centroid; allEccents(frame) = eccent; allMajAxL(frame) = mjaxl; allAspRatio_ml(frame) = aspratio_ml; if do_aspratio_ind allAspRatio_ind(frame) = aspratio_ind; allA0{frame} = A0; allA0c{frame} = A0c; end allLengths(frame) = length(iBdy); allAreas(frame) = area; allBW(:,:,frame) = BW; alliBW(frame,:) = iBW; alloBW(frame,:) = oBW; allinner(frame,:) = inner; allouter(frame,:) = outer; if frame>1 deltaArea = 100*(allAreas(frame,:)/allAreas(frame-1,:)-1); deltaVelocity = norm(allCentroids(frame,:)-allCentroids(frame-1,:)); fprintf('Change area: %6.1f%%, Centroid displacement %6.2f pixels\n',deltaArea,deltaVelocity); else fprintf('\n'); end end % end of nframes %% make kymographs and save data %---------------------------------------------------- dataRaw = v2struct; % pack all variables into single structure array if makekymo kf = kymofixed; % e.g. [2 , 0] where 2 means fixedtheta, and 0 means varying for i = 1:length(kf) binData = cell(1,4); if i==1 [allkymos, alignall, binData1,binData2] = makekymos(dataRaw,kf(i)); else [allkymos, alignall] = makekymos(dataRaw,kf(i)); end tmpKymo.alignangles = binData1(:,1); tmpKymo.alignDtheta = binData1(:,2); tmpKymo.alignbinidx = binData1(:,3); tmpKymo.edges = binData2; % these are the kymos with aligned quantities centered for display tmpKymo.kymovel = allkymos{1}; tmpKymo.kymocurv = allkymos{2}; tmpKymo.kymodcumul = allkymos{3}; tmpKymo.kymofluor = allkymos(4:end); % mch-memb, yfp-actin, cfp-actA % this is never binned tmpKymo.alignbdy = alignall(:,1); % these are the aligned quantities just as lists % these may be spatially binned if using kymofixedtheta or kymofixedpts tmpKymo.alignvel = alignall(:,2); tmpKymo.aligncurv = alignall(:,3); tmpKymo.aligndcumul = alignall(:,4); tmpKymo.alignfluor = alignall(:,5:end); if kf(i) > 0 % kymofixedpts or fixedtheta, vs i == 1 % save kymooutput dataKymoBin = tmpKymo; % plot movies and save only based on the first kymograph set (fixedpts or fixedtheta usually) plottseries(dataParams,dataRaw,makekymo,tmpKymo,kf(i)); % plot figures and output them pause(0.1); figsout0 = plotmovement(dataParams,dataRaw); count = 10; figsoutFix(1) = plotphyskymos(dataParams,tmpKymo,kf(i),count); figsoutFix(2) = plotbiochemkymos(dataParams,tmpKymo,kf(i),count); %figsout0 = v2struct(f3); %figsoutFix = v2struct(f12,f13); else % kymovarying % save kymooutput dataKymoVar = tmpKymo; % plot figures and output them count = 20; figsoutVar(1) = plotphyskymos(dataParams,tmpKymo,kf(i),count); figsoutVar(2) = plotbiochemkymos(dataParams,tmpKymo,kf(i),count); %figsoutVar = v2struct(f22,f23); end end if savefigs == 2 save(strcat(outputfolder,guvname,'_',tstamp1,'_data'),'dataRaw','dataParams',... 'dataKymoBin','dataKymoVar') savefig(figsout0,strcat(outputfolder,guvname,'_',tstamp1,'_figsout0.fig'),'compact'); savefig(figsoutFix,strcat(outputfolder,guvname,'_',tstamp1,'_figsoutFix.fig'),'compact'); savefig(figsoutVar,strcat(outputfolder,guvname,'_',tstamp1,'_figsoutVar.fig'),'compact'); elseif savefigs == 1 save(strcat(outputfolder,guvname,'_',tstamp1,'_data'),'dataParams',... 'dataKymoBin','dataKymoVar'); savefig(figsout0,strcat(outputfolder,guvname,'_',tstamp1,'_figsout0.fig'),'compact'); savefig(figsoutFix,strcat(outputfolder,guvname,'_',tstamp1,'_figsoutFix.fig'),'compact'); savefig(figsoutVar,strcat(outputfolder,guvname,'_',tstamp1,'_figsoutVar.fig'),'compact'); end else % plot tseries that are not based on kymos (like bounding ellipse and alignment and masks) plottseries(dataParams,dataRaw,makekymo); figsout0 = plotmovement(dataParams,dataRaw); if savefigs == 2 save(strcat(outputfolder,guvname,'_',tstamp1,'_data'),'dataRaw','dataParams'); savefig(figsout0,strcat(outputfolder,guvname,'_',tstamp1,'_figsout0.fig'),'compact'); elseif savefigs == 1 save(strcat(outputfolder,guvname,'_',tstamp1,'_data'),'dataParams'); savefig(figsout0,strcat(outputfolder,guvname,'_',tstamp1,'_figsout0.fig'),'compact'); end end end % end of runall %% ------------ auxiliary functions function plottseries(dataParams,dataRaw,makekymo,varargin) if ~isempty(varargin) tmpKymo = varargin{1}; kf = varargin{2}; end cols = dataParams.cols; rows = dataParams.rows; nframes = dataParams.nframes; listFrames = dataParams.listFrames; fnames = dataParams.fnames; % name of channels membch = dataParams.membch; fluorch = dataParams.fluorch; guvname = dataParams.guvname; savgolay_smooth = dataParams.savgolay_smooth; computecurv = dataParams.computecurv; im0 = dataParams.im0; imL = dataParams.imL; do_plotEllipseries = dataParams.do_plotEllipseries; do_plotTimeseries = dataParams.do_plotTimeseries; do_plotAlignseries = dataParams.do_plotAlignseries; do_plotMaskseries = dataParams.do_plotMaskseries; do_aspratio_ind = dataParams.do_aspratio_ind; savefigs = dataParams.savefigs; tstamp1 = dataParams.tstamp1; outputfolder = dataParams.outputfolder; sat = dataParams.sat; satlow = dataParams.satlow; deltaf = dataParams.deltaf; nfluor = dataParams.nfluor; tickstep = dataParams.tickstep; minfont = dataParams.minfont; basystem = dataParams.basystem; if makekymo % these are the kymos kymovel = tmpKymo.kymovel; % aligned velocities kymodcumul = tmpKymo.kymodcumul; kymocurv = tmpKymo.kymocurv; % aligned curvatures kymofluor = tmpKymo.kymofluor; % aligned kymos for fluorescent channels % aligned data, these can be spatially binned, but not centered for display like kymos alignbdy = tmpKymo.alignbdy; % this will never be binned - longer than alignvel/curv/fluor alignvel = tmpKymo.alignvel; aligndcumul = tmpKymo.aligndcumul; aligncurv = tmpKymo.aligncurv; alignfluor = tmpKymo.alignfluor; % % full aligned data, with varying # points (m1*perim) or fixed long perimeter (max m1*perim) % longbdy = dataKymos.longbdy; % longvel = dataKymos.longvel; % longcurv = dataKymos.longcurv; % longdcumul = dataKymos.longdcumul; % longfluor = dataKymos.longfluor; end cent = dataRaw.allCentroids; areas = dataRaw.allAreas; allImages = dataRaw.allImages; allIrgb = dataRaw.allIrgb; alloBW = dataRaw.alloBW; alliBW = dataRaw.alliBW; allinner = dataRaw.allinner; allouter = dataRaw.allouter; if do_aspratio_ind % these are aspect ratio quantities allA0 = dataRaw.allA0; allA0c = dataRaw.allA0c; end allBoundaries = dataRaw.allBoundaries; %vflag = any(~isnan(kymovel(:))); % check if velocity is available %% plot min bounding ellipse if do_plotEllipseries if ~do_aspratio_ind % must have do_aspratio_ind OFF to compute aspect ratio quantities disp('No bounding ellipse available: turn do_aspratio_ind ON for independent aspect ratios with A0 and A0c') else pause(0.1) % fname = strcat(outputfolder,guvname,'_',tstamp1,'_bounding_ellipse'); % if strcmp('linux',basystem) % v = VideoWriter(fname,'Motion JPEG AVI'); % else % v = VideoWriter(fname,'MPEG-4'); % end % % v.FrameRate = 5; %playback framerate % open(v); for j = 1:nframes movieFrame=listFrames(j); %im0:imL; % TIFF indices tmpIrgb = allIrgb{j,membch}; bdyraw = allBoundaries{j}; tmpA0 = allA0{j}; tmpA0c = allA0c{j}; byraw = bdyraw(:,1); bxraw = bdyraw(:,2); xpts = [bxraw,byraw]; f24 = figure(24); imshow(tmpIrgb); axis equal; hold on; plot(xpts(:,1),xpts(:,2),'r','LineWidth',3); FontSize = minfont; %min(minfont,floor(min(cols,rows)/10)); text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) %set(gca,'YDir','reverse'); Ellipse_plot(tmpA0,tmpA0c); hold off; drawnow % save as tiff series if savefigs > 0 pause(0.05) F24 = getframe(f24); if movieFrame == im0 imwrite(F24.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_bounding_ellipse.tif')); else imwrite(F24.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_bounding_ellipse.tif'),... 'WriteMode','append'); end end %writeVideo(v,F4.cdata); end %close(v); close (f24); end end %% plot boundary alignment if do_plotAlignseries for j = 1:nframes % indices for arrays movieFrame=listFrames(j); %im0:imL; % TIFF indices f60 = figure(60); tmpbdy = alignbdy{j}; tmpI1rgb = allIrgb{j,membch}; x0 = cent(j,1); y0 = cent(j,2); imshow(tmpI1rgb); axis equal; hold on; ax60 = plot(x0,y0,'m*','MarkerSize',12,'LineWidth',3); ax60.Parent.Title.String = strcat(guvname,{' '},'tracking reference angle (center of kymograph)'); FontSize = minfont; %min(minfont,floor(min(cols,rows)/10)); text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) %plot(tmpbdy(:,2),tmpbdy(:,1),'ro','MarkerSize',10); plot(tmpbdy(:,2),tmpbdy(:,1),'r-','LineWidth',3); midpoint = round(length(tmpbdy)/2); plot(tmpbdy(midpoint,2),tmpbdy(midpoint,1),'c*','MarkerSize',12,'LineWidth',3); hold off; drawnow % save as tiff series if savefigs > 0 pause(0.05) F60 = getframe(f60); if movieFrame == im0 imwrite(F60.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_align.tif')); else imwrite(F60.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_align.tif'),... 'WriteMode','append'); end end end close (f60); end % end of plot alignseries %% timeseries tiff main boundary measurements if do_plotTimeseries && makekymo pause(1); % find caxis lims for fluor quants fluorcax = cell(1,nfluor); for k = 1:nfluor if kf > 0 %kymofixedpts || kymofixedtheta a3 = cat(2,alignfluor{:,k}); % these are same dimensions so can concatenize else a3 = cat(2,kymofluor{k}); % can only concat if same dimensions, so use kymo which has padding end fluormax = sat*max(a3(:)); fluormin = satlow*min(a3(:)); fluorcax{k} = [fluormin, fluormax]; end % find caxis lims for bdy quants avel = cat(2,alignvel{:}); adc = cat(2,aligndcumul{:}); acurv = cat(2,aligncurv{:}); velmax = max(abs(avel(:))); dcmax = max(abs(adc(:))); curvmax = max(abs(acurv(:))); velcax = [-velmax,velmax]; dccax = [-dcmax,dcmax]; curvcax = [-curvmax,curvmax]; %fluor quantitites panel flnames = {'membrane marker';'actin';'ActA'}; phnames = {'velocity';'cuvature';'cumulative disp.'}; figpan = struct('names',{flnames(:),phnames(:)},... 'vals',{{alignfluor},{alignvel;aligndcumul;aligncurv}},... 'caxes',{fluorcax(:),{velcax;dccax;curvcax}}); % now plot multi-panel time series together for k = 1:length(figpan) for j = 1:nframes % indices for arrays movieFrame=listFrames(j); %im0:imL; % TIFF indices tmpI1rgb = allIrgb{j,membch}; tmpbdy = alignbdy{j}; tmpvel = alignvel{j}; tmpdc = aligndcumul{j}; tmpcurv = aligncurv{j}; tmpfluor = alignfluor{j,k}; tmpI2rgb = allIrgb{j,k}; fignum = 2; f2vars = v2struct(savgolay_smooth,computecurv,movieFrame,panelnames,panelcaxes,tmpbdy); f2 = displaymultipanel(fignum,f2vars,tmpI1rgb,tmpI1pan); % if k == 2 && movieFrame==im0 % only for actin show triple time series % % plot % f2vars = v2struct(savgolay_smooth,computecurv,movieFrame,panelnames,minfont,velcax,curvcax,panelcaxes); % f2 = displaymultipanel(tmpI1rgb,tmpI2rgb,f2vars,tmpbdy,tmpcurv,tmpfluor); % elseif k == 2 && movieFrame > im0 % velocities available % % then plot % f2vars = v2struct(savgolay_smooth,computecurv,movieFrame,panelnames,minfont,velcax,curvcax,panelcaxes); % f2 = displaymultipanel(tmpI1rgb,tmpI2rgb,f2vars,tmpbdy,tmpcurv,tmpfluor,tmpvel); % else % for other fluor channels, just show the intensity around boundary % f2vars = v2struct(savgolay_smooth,computecurv,movieFrame,panelnames,minfont,velcax,curvcax,panelcaxes); % f2 = displaymultipanel(tmpI1rgb,tmpI2rgb,f2vars,tmpbdy,tmpfluor); % end % if savefigs > 0 % create tiff series %pause(0.05) F2 = getframe(f2); if movieFrame == im0 imwrite(F2.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_main_',fnames{k}(5:end),'.tif')); else imwrite(F2.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_main_',fnames{k}(5:end),'.tif'),... 'WriteMode','append'); end end end % end of nframes %close(f2) end % end of nfluor end %end of plot main time series %% plot mask time series if do_plotMaskseries for k = 1:nfluor % include membrane channel here f8 = figure(8); clf(f8); f8.Position(3:4) = [400, 400]; tiledlayout(2,2,'TileSpacing','compact','Padding','compact'); for j = 1:nframes movieFrame = listFrames(j); nexttile(1,[2,2]); Iflrgb = allIrgb{j,k}; iBW = alliBW{j,k}; oBW = alloBW{j,k}; inner = allinner{j,k}; outer = allouter{j,k}; if makekymo tmpbdy = alignbdy{j}; else tmpbdy = allBoundaries{j}; end f8ax = imshow(Iflrgb); axis equal; hold on; green = cat(3, zeros(rows,cols), ones(rows,cols), zeros(rows,cols)); h8 = imshow(green); % Use mask as the AlphaData for the solid green image. set(h8, 'AlphaData', 0.5*oBW) ; red = cat(3, ones(rows,cols), zeros(rows,cols), zeros(rows,cols)); h9 = imshow(red); % Use mask as the AlphaData for the solid green image. set(h9, 'AlphaData', 0.5*iBW) ; % plot boundary points h10=plot(inner(:,1),inner(:,2),'r--'); plot(outer(:,1),outer(:,2),'c--'); plot(tmpbdy(:,2),tmpbdy(:,1),'w--'); hold off; FontSize = minfont; %min(minfont,floor(min(cols,rows)/10)); text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) h10.Parent.YLim=[1 rows]; h10.Parent.XLim=[1 cols]; title(['Normalized',' ',fnames{k}(5:end)],'Interpreter','none') % f8.Position(3:4) = [400, 400]; drawnow if savefigs>0 %pause(0.05) F8 = getframe(f8); if movieFrame == im0 imwrite(F8.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_mask_',fnames{k}(5:end),'.tif')); else imwrite(F8.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_mask_',fnames{k}(5:end),'.tif'),... 'WriteMode','append'); end end end end end %if do_plotmaskseries end function figsout = plotphyskymos(dataParams,tmpKymo,kf,count) cols = dataParams.cols; rows = dataParams.rows; nframes = dataParams.nframes; listFrames = dataParams.listFrames; fnames = dataParams.fnames; % name of channels membch = dataParams.membch; fluorch = dataParams.fluorch; guvname = dataParams.guvname; savgolay_smooth = dataParams.savgolay_smooth; computecurv = dataParams.computecurv; im0 = dataParams.im0; imL = dataParams.imL; savefigs = dataParams.savefigs; tstamp1 = dataParams.tstamp1; outputfolder = dataParams.outputfolder; t_rapa = dataParams.t_rapa; tickstep = dataParams.tickstep; dt = dataParams.dt; deltaf = dataParams.deltaf; nbins = dataParams.nbins; kymovel = tmpKymo.kymovel; % aligned velocities kymocurv = tmpKymo.kymocurv; % aligned curvatures kymodcumul = tmpKymo.kymodcumul; % aligned cumulative displacements kymofluor = tmpKymo.kymofluor; % aligned kymos for fluorescent channels alignbdy = tmpKymo.alignbdy; %allBoundaries; alignvel = tmpKymo.alignvel; %allVelocities; aligncurv = tmpKymo.aligncurv; alignfluor = tmpKymo.alignfluor; %allIntensities; edges = tmpKymo.edges; % % cent = dataRaw.allCentroids; % areas = dataRaw.allAreas; % allImages = dataRaw.allImages; % allIrgb = dataRaw.allIrgb; vflag = any(~isnan(kymovel(:))); % check if velocity is available tfs = 17; %title font size afs = 14; % axes labels font size xlinec = [0.1 0.1 0.1]; %% ---------------- plot kymographs f12 = figure(2+count); t2 = tiledlayout('flow'); %t2=tiledlayout(3,3); t2.Padding = 'compact'; t2.TileSpacing = 'compact'; nexttile(1,[1 3]) if vflag ax4 = imagesc(kymovel,'AlphaData',~isnan(kymovel)); cmax = max(abs(kymovel(:))); colormap(ax4.Parent,'parula'); ax4.Parent.CLim = [-cmax,cmax]; colorbar; % xticks(1:tickstep:size(kymovel,2)); xticklabels({im0:tickstep:imL}); xticks(tickstep:tickstep:size(kymovel,2)+1); xticklabels({im0-1+tickstep:tickstep:imL}); % recall: % edges = - pi : 2*pi/nbins : pi ; % edges = 1 : (n-1)/nbins : n; hold on xline(t_rapa+0.5,'LineWidth',3,'Color',xlinec,'LineStyle','--'); if guvname(1) == 'L' if kf == 2 %fixedtheta plot(im0-0.5,0,'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); elseif kf == 1 %fixedpoints plot(im0-0.5,round(nbins/2),'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); end end hold off if kf == 2 %kymofixedtheta edgeticks = 0.5:nbins/6:nbins+0.5; % edgelabs = pi : -2*pi/6 : -pi; edgelabs = {'\pi';'2\pi/3';'\pi/3';'0';'-\pi/3';'-2\pi/3';'-\pi'}; yticks(edgeticks); yticklabels(edgelabs); % elseif kymofixedpts % edgeticks = 0.5:nbins/6:nbins+0.5; % edgelabs = 0:nbins/6:nbins; % yticks(edgeticks); yticklabels({edgelabs}); end ax4.Parent.XLabel.String = 'Frame Number'; %ax4.Parent.XLabel.FontSize = afs; if kf == 2 %kymofixedtheta ax4.Parent.YLabel.String = 'Angular position (binned)'; elseif kf == 1 %kymofixedpts ax4.Parent.YLabel.String = 'Perimeter (binned points)'; else ax4.Parent.YLabel.String = 'Perimeter (points)'; end ax4.Parent.YLabel.FontName = 'Arial'; %ax4.Parent.YLabel.FontSize = afs; shading flat; end ax4.Parent.FontSize = afs; title('Velocity (pixels/frame) along perimeter','FontSize',tfs); %------------------------------------ nexttile(4,[1 3]) ax5 = imagesc(kymodcumul,'AlphaData',~isnan(kymodcumul)); shading flat % cm = zeros(256,3,'uint8'); % cm(1:256,2) = linspace(0,255,256); % cm(1:256,3) = linspace(0,255,256); cmax = max(abs(kymodcumul(:))); %cmin = min(kymo2(kymo2>0)); colormap(ax5.Parent,'parula'); %ax5.Parent.Colormap = cm; ax5.Parent.CLim = [-cmax,cmax]; colorbar; % xticks(1:tickstep:size(kymofluor,2)); xticklabels({im0:tickstep:imL}); xticks(tickstep:tickstep:size(kymodcumul,2)); xticklabels({im0-1+tickstep:tickstep:imL}); hold on xline(t_rapa+0.5,'LineWidth',3,'Color',xlinec,'LineStyle','--'); if guvname(1) == 'L' if kf == 2 %fixedtheta plot(im0-0.5,0,'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); elseif kf == 1 %fixedpoints plot(im0-0.5,round(nbins/2),'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); end end hold off if kf == 2 %kymofixedtheta edgeticks = 0.5:nbins/6:nbins+0.5; % edgelabs = pi : -2*pi/6 : -pi; edgelabs = {'\pi';'2\pi/3';'\pi/3';'0';'-\pi/3';'-2\pi/3';'-\pi'}; yticks(edgeticks); yticklabels(edgelabs); end ax5.Parent.XLabel.String = 'Frame Number'; ax5.Parent.XLabel.FontName = 'Arial'; %ax5.Parent.XLabel.FontSize = afs; if kf ==2 %kymofixedtheta ax5.Parent.YLabel.String = 'Angular position (binned)'; elseif kf == 1 %kymofixedpts ax5.Parent.YLabel.String = 'Perimeter (binned points)'; else ax5.Parent.YLabel.String = 'Perimeter (points)'; end %ax5.Parent.YLabel.FontName = 'Arial'; ax5.Parent.YLabel.FontSize = afs; ax5.Parent.FontSize = afs; title('Cumulative displacement (pixels) along perimeter','FontSize',tfs); %---------------------------------------------- nexttile(7,[1 3]) ax6 = imagesc(kymocurv,'AlphaData',~isnan(kymocurv)); shading flat % cmax = max(abs(kymocurv(:))); % colormap(ax6.Parent,'parula'); % ax6.Parent.CLim = [-cmax,cmax]; % colorbar; % BINARIZATION % clrs = [0.2422 0.1504 0.6603; 0.9290 0.6940 0.1250]; % ax6.Parent.CLim = [0,1]; % colormap(ax6.Parent,clrs); % colorbar; % ORIGINAL COLOR SCHEME cmax = max(abs(kymocurv(:))); colormap(ax6.Parent,'parula'); ax6.Parent.CLim = [-cmax,cmax]; colorbar; % xticks(1:tickstep:size(kymocurv,2)); xticklabels({im0:tickstep:imL}); xticks(tickstep:tickstep:size(kymocurv,2)); xticklabels({im0-1+tickstep:tickstep:imL}); hold on xline(t_rapa+0.5,'LineWidth',3,'Color',xlinec,'LineStyle','--'); if guvname(1) == 'L' if kf == 2 %fixedtheta plot(im0-0.5,0,'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); elseif kf == 1 %fixedpoints plot(im0-0.5,round(nbins/2),'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); end end hold off if kf == 2 %kymofixedtheta edgeticks = 0.5:nbins/6:nbins+0.5; %edgelabs = pi : -2*pi/6 : -pi; edgelabs = {'\pi';'2\pi/3';'\pi/3';'0';'-\pi/3';'-2\pi/3';'-\pi'}; yticks(edgeticks); yticklabels(edgelabs); end ax6.Parent.XLabel.String = 'Frame Number'; ax6.Parent.XLabel.FontName = 'Arial'; %ax6.Parent.XLabel.FontSize = afs; if kf == 2 %kymofixedtheta ax6.Parent.YLabel.String = 'Angular position (binned)'; elseif kf == 1 %kymofixedpts ax6.Parent.YLabel.String = 'Perimeter (binned points)'; else ax6.Parent.YLabel.String = 'Perimeter (points)'; end ax6.Parent.YLabel.FontName = 'Arial'; %ax6.Parent.YLabel.FontSize = afs; ax6.Parent.FontSize = afs; if computecurv == 1 title('Curvature (circumcenter) along perimeter','FontSize',tfs); else title('Curvature (fitcircle) along perimeter','FontSize',tfs); end % title(t2,[guvname,': physical kymographs']); % adjust the figure positions %[left bottom width height] f12.Position(3:4) = [700 700]; drawnow; % %% adjust font size in all figures % fh = findall(0,'Type','Figure'); % set( findall(fh, '-property', 'fontsize'), 'fontsize', 14) %figsout = v2struct(f12); figsout = f12; if savefigs>0 && kf>0 print(f12,strcat(outputfolder,guvname,'_',tstamp1,'_fig_physkymos_Fix'),'-dpng') elseif savefigs>0 && kf==0 print(f12,strcat(outputfolder,guvname,'_',tstamp1,'_fig_physkymos_Var'),'-dpng') end end function figsout = plotbiochemkymos(dataParams,tmpKymo,kf,count) cols = dataParams.cols; rows = dataParams.rows; nframes = dataParams.nframes; listFrames = dataParams.listFrames; fnames = dataParams.fnames; % name of channels membch = dataParams.membch; fluorch = dataParams.fluorch; guvname = dataParams.guvname; savgolay_smooth = dataParams.savgolay_smooth; computecurv = dataParams.computecurv; im0 = dataParams.im0; imL = dataParams.imL; savefigs = dataParams.savefigs; tstamp1 = dataParams.tstamp1; outputfolder = dataParams.outputfolder; t_rapa = dataParams.t_rapa; tickstep = dataParams.tickstep; dt = dataParams.dt; deltaf = dataParams.deltaf; nbins = dataParams.nbins; sat = dataParams.sat; satlow = dataParams.satlow; kymovel = tmpKymo.kymovel; % aligned velocities kymocurv = tmpKymo.kymocurv; % aligned curvatures kymodcumul = tmpKymo.kymodcumul; % aligned cumulative displacements kymofluor = tmpKymo.kymofluor; % aligned kymos for fluorescent channels alignbdy = tmpKymo.alignbdy; %allBoundaries; alignvel = tmpKymo.alignvel; %allVelocities; aligncurv = tmpKymo.aligncurv; alignfluor = tmpKymo.alignfluor; %allIntensities; edges = tmpKymo.edges; vflag = any(~isnan(kymovel(:))); % check if velocity is available tfs = 17; % title font size afs = 14; % axes font size xlinec = [0.1 0.1 0.1]; %% ---------------- plot kymographs f13 = figure(3+count); t2 = tiledlayout('flow'); %t2=tiledlayout(3,3); t2.Padding = 'compact'; t2.TileSpacing = 'compact'; nexttile(1,[1 3]) if vflag kymfl = kymofluor{fluorch(1)}; ax4 = imagesc(kymfl,'AlphaData',~isnan(kymfl)); cmax = max(kymfl(:)); cmin = min(kymfl(:)); colormap(ax4.Parent,'parula'); ax4.Parent.CLim = [satlow*cmin,sat*cmax]; colorbar; % xticks(1:tickstep:size(kymovel,2)); xticklabels({im0:tickstep:imL}); xticks(tickstep:tickstep:size(kymfl,2)+1); xticklabels({im0-1+tickstep:tickstep:imL}); hold on xline(t_rapa+0.5,'LineWidth',3,'Color',xlinec,'LineStyle','--'); if guvname(1) == 'L' if kf == 2 %fixedtheta plot(im0-0.5,0,'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); elseif kf == 1 %fixedpoints plot(im0-0.5,round(nbins/2),'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); end end hold off if kf == 2 %kymofixedtheta edgeticks = 0.5:nbins/6:nbins+0.5; % edgelabs = pi : -2*pi/6 : -pi; edgelabs = {'\pi';'2\pi/3';'\pi/3';'0';'-\pi/3';'-2\pi/3';'-\pi'}; yticks(edgeticks); yticklabels(edgelabs); end ax4.Parent.XLabel.String = 'Frame Number'; %ax4.Parent.XLabel.FontSize = afs; if kf == 2 %kymofixedtheta ax4.Parent.YLabel.String = 'Angular position (binned)'; elseif kf == 1 % kymofixedpts ax4.Parent.YLabel.String = 'Perimeter (binned points)'; else ax4.Parent.YLabel.String = 'Perimeter (points)'; end ax4.Parent.YLabel.FontName = 'Arial'; %ax4.Parent.YLabel.FontSize = afs; shading flat; end ax4.Parent.FontSize = afs; title('Membrane marker intensity along perimeter','FontSize',tfs); % ------------------------ nexttile(4,[1 3]) kymfl = kymofluor{fluorch(3)}; ax6 = imagesc(kymfl,'AlphaData',~isnan(kymfl)); shading flat cmax = max(kymfl(:)); cmin = min(kymfl(:)); colormap(ax6.Parent,'parula'); ax6.Parent.CLim = [satlow*cmin,sat*cmax]; colorbar; % xticks(1:tickstep:size(kymocurv,2)); xticklabels({im0:tickstep:imL}); xticks(tickstep:tickstep:size(kymfl,2)); xticklabels({im0-1+tickstep:tickstep:imL}); hold on xline(t_rapa+0.5,'LineWidth',3,'Color',xlinec,'LineStyle','--'); if guvname(1) == 'L' if kf == 2 %fixedtheta plot(im0-0.5,0,'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); elseif kf == 1 %fixedpoints plot(im0-0.5,round(nbins/2),'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); end end hold off if kf == 2 %kymofixedtheta edgeticks = 0.5:nbins/6:nbins+0.5; % edgelabs = pi : -2*pi/6 : -pi; edgelabs = {'\pi';'2\pi/3';'\pi/3';'0';'-\pi/3';'-2\pi/3';'-\pi'}; yticks(edgeticks); yticklabels(edgelabs); end ax6.Parent.XLabel.String = 'Frame Number'; ax6.Parent.XLabel.FontName = 'Arial'; %ax6.Parent.XLabel.FontSize = afs; if kf == 2 %kymofixedtheta ax6.Parent.YLabel.String = 'Angular position (binned)'; elseif kf == 1 %kymofixedpts ax6.Parent.YLabel.String = 'Perimeter (binned points)'; else ax6.Parent.YLabel.String = 'Perimeter (points)'; end ax6.Parent.YLabel.FontName = 'Arial'; %ax6.Parent.YLabel.FontSize = afs; ax6.Parent.FontSize = afs; title('ActA intensity along perimeter','FontSize',tfs); % %------------------------------ nexttile(7,[1 3]) kymfl = kymofluor{fluorch(2)}; ax5 = imagesc(kymfl,'AlphaData',~isnan(kymfl)); shading flat % cmax = max(kymfl(:)); % cmin = min(kymfl(:)); % colormap(ax5.Parent,'parula'); % ax5.Parent.CLim = [satlow*cmin,sat*cmax]; % colorbar; % BINARIZATION % clrs = [0.2422 0.1504 0.6603; 0.9290 0.6940 0.1250]; % ax5.Parent.CLim = [0,1]; % colormap(ax5.Parent,clrs); % colorbar; % ORIGINAL COLOR SCHEME cmax = max(kymfl(:)); cmin = min(kymfl(:)); colormap(ax5.Parent,'parula'); ax5.Parent.CLim = [satlow*cmin,sat*cmax]; colorbar; % xticks(1:tickstep:size(kymofluor,2)); xticklabels({im0:tickstep:imL}); xticks(tickstep:tickstep:size(kymfl,2)); xticklabels({im0-1+tickstep:tickstep:imL}); hold on xline(t_rapa+0.5,'LineWidth',3,'Color',xlinec,'LineStyle','--'); if guvname(1) == 'L' if kf == 2 %fixedtheta plot(im0-0.5,0,'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); elseif kf == 1 %fixedpoints plot(im0-0.5,round(nbins/2),'r>','MarkerSize',8,'MarkerFaceColor',[1 0.5 0]); end end hold off if kf == 2 %kymofixedtheta edgeticks = 0.5:nbins/6:nbins+0.5; %edgelabs = pi : -2*pi/6 : -pi; edgelabs = {'\pi';'2\pi/3';'\pi/3';'0';'-\pi/3';'-2\pi/3';'-\pi'}; yticks(edgeticks); yticklabels(edgelabs); end ax5.Parent.XLabel.String = 'Frame Number'; ax5.Parent.XLabel.FontName = 'Arial'; %ax5.Parent.XLabel.FontSize = afs; if kf == 2 %kymofixedtheta ax5.Parent.YLabel.String = 'Angular position (binned)'; elseif kf == 1 %kymofixedpts ax5.Parent.YLabel.String = 'Perimeter (binned points)'; else ax5.Parent.YLabel.String = 'Perimeter (points)'; end ax5.Parent.YLabel.FontName = 'Arial'; %ax5.Parent.YLabel.FontSize = afs; ax5.Parent.FontSize = afs; title('Actin intensity along perimeter','FontSize',tfs); % % title(t2,[guvname,': biochemical kymographs (top ',num2str(100*(1-sat)),... % '% sat.)']); title(t2,[guvname,': biochemical kymographs']); drawnow % %% adjust font size in all figures % fh = findall(0,'Type','Figure'); % set( findall(fh, '-property', 'fontsize'), 'fontsize', 14) % adjust the figure positions %[left bottom width height] f13.Position(3:4) = [700 700]; %figsout = v2struct(f13); figsout = f13; if savefigs>0 && kf>0 print(f13,strcat(outputfolder,guvname,'_',tstamp1,'_fig_biochemkymos_Fix'),'-dpng') elseif savefigs>0 && kf==0 print(f13,strcat(outputfolder,guvname,'_',tstamp1,'_fig_biochemkymos_Var'),'-dpng') end end function figsout = plotmovement(dataParams,dataRaw) cols = dataParams.cols; rows = dataParams.rows; nframes = dataParams.nframes; listFrames = dataParams.listFrames; guvname = dataParams.guvname; im0 = dataParams.im0; imL = dataParams.imL; savefigs = dataParams.savefigs; tstamp1 = dataParams.tstamp1; outputfolder = dataParams.outputfolder; tickstep = dataParams.tickstep; t_rapa = dataParams.t_rapa; do_aspratio_ind = dataParams.do_aspratio_ind; do_fixedbox = dataParams.do_fixedbox; do_plotBoundaryOverlay = dataParams.do_plotBoundaryOverlay; cent = dataRaw.allCentroids; areas = dataRaw.allAreas; aspratio_ml = dataRaw.allAspRatio_ml; aspratio_ind = dataRaw.allAspRatio_ind; eccent = dataRaw.allEccents; bdy = dataRaw.allBoundaries; allMajAxL = dataRaw.allMajAxL; afs = 14; % axes font size tfs = afs + 3; % title font size %% ------------- plot centroid and velocity f3 = figure(3); t1=tiledlayout(4,3); %(5,3) if including centroid velocity fdelta = imL-im0+1; cm=colormap(parula(fdelta)); % use as many colors as there frames % minRow = rows; % maxRow = 1; % minCol = cols; % maxCol = 1; minRow = 1; maxRow = rows; minCol = 1; maxCol = cols; velocityCentroid = [nan;vecnorm(diff(cent)')']; if do_fixedbox % FOR FIXED BOX SIZE boxt1 = [1,180,1,180]; % xmin, xmax, ymin, ymax maxrc = boxt1([2,4]); dispadj = round(maxrc/2); % center x, center y dxadj = cent(1,:)-fliplr(dispadj); % x, y adjust else boxt1 = [minCol maxCol minRow maxRow]; dxadj = [0 0]; end for k = 1:nframes % indices for arrays movieFrame=listFrames(k); %movieframe, usually same as im0:imL; % TIFF indices %b = mod(a,m) returns the remainder after division of a by m, %thiscm = cm(mod(k-1,255)+1,:); thiscm = cm(movieFrame-im0+1,:); thiscm2 = [0 0.4470 0.7410]; backcm = 0.3*[1 1 1]; tmpbdy=bdy{k}; if do_fixedbox boxsize = boxt1; else % adjust overall min, max after each frame minCol = min(minCol,min(tmpbdy(:,2)-dxadj(1))); maxCol = max(maxCol,max(tmpbdy(:,2)-dxadj(1))); minRow = min(minRow,min(tmpbdy(:,1)-dxadj(2))); maxRow = max(maxRow,max(tmpbdy(:,1)-dxadj(2))); boxsize = [minCol maxCol minRow maxRow]; end nexttile(1) % plot centroid: already in x, y coords (not row, col coords) %ax1=plot(cent(k,1),rows+1-cent(k,2),'.','Color',thiscm,'MarkerSize',30); ax1=plot(cent(k,1)-dxadj(1),cent(k,2)-dxadj(2),'.','Color',thiscm,'MarkerSize',30); axis equal, hold on set(ax1.Parent,'Color',backcm) set(ax1.Parent,'Ydir','reverse') axis(boxsize) nexttile(2) % plot boundary: convert row index to y coord %ax2=plot(tmpbdy(:,2),rows+1-tmpbdy(:,1),'-','Color',thiscm,'LineWidth',2); ax2=plot(tmpbdy(:,2)-dxadj(1),tmpbdy(:,1)-dxadj(2),'-','Color',thiscm,'LineWidth',2); axis equal, hold on set(ax2.Parent,'Color',backcm) set(ax2.Parent,'Ydir','reverse') axis(boxsize) ax2.Parent.CLim = [im0 imL+1]; end % add colorbar to ax2 %cc = colorbar('Ticks',im0+0.5:tickstep:imL+0.5,'TickLabels',{im0:tickstep:imL}); cc = colorbar('Ticks',im0-0.5+tickstep:tickstep:imL+0.5,'TickLabels',{im0-1+tickstep:tickstep:imL}); cc.Label.String = 'Frames'; cc.Label.FontSize = afs; cc.FontSize = afs; % next plots are stair plots with only 1 color, FAST % stair eccentricity, major axis length/asp ratio, area, centroid velocity plot with color variation nexttile(4,[1 3]) % plot eccentricity tmpeccent = nan(fdelta,1); tmpeccent(listFrames)=eccent; ax3a = stairs([tmpeccent;tmpeccent(end)],'Color',thiscm2,'LineWidth',2); ax3a.Parent.CLim = [im0 imL+1]; %set(gca,'Color',backcm) hold on nexttile(7,[1 3]) % % plot major axis length % tmpmaxl = nan(fdelta,1); % tmpmaxl(listFrames)=allMajAxL; % ax3b=stairs([tmpmaxl;tmpmaxl(end)],'.-','Color',thiscm2,'LineWidth',2); % ax3b.Parent.CLim = [im0 imL+1]; % hold on % plot aspect ratio tmpaspr = nan(fdelta,1); % if do_aspratio_ind % tmpaspr(listFrames)=aspratio_ind; % else % use matlab regionprops % tmpaspr(listFrames)=aspratio_ml; % end tmpaspr(listFrames)=aspratio_ml; ax3b = stairs([tmpaspr;tmpaspr(end)],'Color',thiscm2,'LineWidth',2); ax3b.Parent.CLim = [im0 imL+1]; %set(gca,'Color',backcm) hold on nexttile(10,[1 3]) % plot areas tmpareas = nan(fdelta,1); tmpareas(listFrames)=areas; ax3l=stairs([tmpareas;tmpareas(end)],'.-','Color',thiscm2,'LineWidth',2); ax3l.Parent.CLim = [im0 imL+1]; %set(gca,'Color',backcm) hold on % nexttile(13,[1 3]) % % stair centroid velocity plot with color variation % tmpcvel = nan(fdelta,1); % tmpcvel(listFrames)=velocityCentroid; % ax3r=stairs([tmpcvel;tmpcvel(end)],'.-','Color',thiscm2,'LineWidth',2); % ax3r.Parent.CLim = [im0 imL+1]; % %set(gca,'Color',backcm) % hold on if do_fixedbox boxsize = boxt1; else %extra padding when adjusting box for every GUV minCol = max( 1,minCol-5); maxCol = min(cols,maxCol+5); minRow = max( 1,minRow-5); maxRow = min(rows,maxRow+5); boxsize = [minCol maxCol minRow maxRow]; end % nexttile(1) axis(boxsize), axis equal hold off ax1.Parent.FontName = 'Arial'; ax1.Parent.FontSize = afs; ax1.Parent.XLim = boxsize(1:2); ax1.Parent.YLim = boxsize(3:4); title('Centroid Position','FontSize',tfs) %ax1.Parent.YTickLabel = ax1.Parent.YTickLabel(end:-1:1); % nexttile(2) axis(boxsize), axis equal hold off ax2.Parent.FontName = 'Arial'; ax2.Parent.FontSize = afs; ax2.Parent.XLim = boxsize(1:2); ax2.Parent.YLim = boxsize(3:4); title('Boundary','FontSize',tfs) %ax2.Parent.YTickLabel = ax2.Parent.YTickLabel(end:-1:1); % nexttile(4,[1 3]) xline(t_rapa+1,'LineWidth',2) hold off ax3a.Parent.YLabel.String = 'Eccentricity'; ax3a.Parent.YLim=[0 1]; ax3a.Parent.YLabel.FontName = 'Arial'; %ax3a.Parent.YLabel.FontSize = afs; ax3a.Parent.XLim=[im0 imL+1]; ax3a.Parent.XTick=im0-0.5+tickstep:tickstep:imL+0.5; ax3a.Parent.XTickLabel=im0-1+tickstep:tickstep:imL; ax3a.Parent.XLabel.String='Frames'; %ax3a.Parent.XLabel.FontSize=afs; ax3a.Parent.FontSize = afs; % nexttile(7,[1 3]) xline(t_rapa+1,'LineWidth',2) hold off ax3b.Parent.YLabel.String = 'Aspect Ratio'; ax3b.Parent.YLim(1) = 1; % lower limit of aspect ratios %ax3b.Parent.YLim=[1 3]; % range of aspect ratios %ax3b.Parent.YLabel.String = 'Major Axis Length (pixels)'; ax3b.Parent.YLabel.FontName = 'Arial'; %ax3b.Parent.YLabel.FontSize = afs; ax3b.Parent.XLim=[im0 imL+1]; ax3b.Parent.XTick=im0-0.5+tickstep:tickstep:imL+0.5; ax3b.Parent.XTickLabel=im0-1+tickstep:tickstep:imL; ax3b.Parent.XLabel.String='Frames'; %ax3b.Parent.XLabel.FontSize=afs; ax3b.Parent.FontSize = afs; % nexttile(10,[1 3]) xline(t_rapa+1,'LineWidth',2) hold off ax3l.Parent.YLabel.String = 'Area (pixels)'; ax3l.Parent.YLabel.FontName = 'Arial'; %ax3l.Parent.YLabel.FontSize = afs; ax3l.Parent.XLim=[im0 imL+1]; ax3l.Parent.XTick=im0-0.5+tickstep:tickstep:imL+0.5; ax3l.Parent.XTickLabel=im0-1+tickstep:tickstep:imL; ax3l.Parent.XLabel.String='Frames'; %ax3l.Parent.XLabel.FontSize=afs; ax3l.Parent.FontSize = afs; % % nexttile(13,[1 3]) % xline(t_rapa+1,'LineWidth',2) % hold off % ax3r.Parent.YLabel.String = 'Centroid Velocity (pixels/frame)'; % ax3r.Parent.YLabel.FontName = 'Arial'; % %ax3r.Parent.YLabel.FontSize = afs; % ax3r.Parent.XLim=[im0 imL+1]; % ax3r.Parent.XTick=im0-0.5+tickstep:tickstep:imL+0.5; % ax3r.Parent.XTickLabel=im0-1+tickstep:tickstep:imL; % ax3r.Parent.XLabel.String='Frames'; % %ax3r.Parent.XLabel.FontSize=afs; % ax3r.Parent.FontSize = afs; % t1.Padding = 'compact'; t1.TileSpacing = 'compact'; t1.Title.String = guvname; t1.Title.FontSize = 12; t1.Title.FontName = 'Arial'; t1.Title.Color = [0 0 0]; % adjust the figure positions %[left bottom width height] f3.Position(3:4) = [800, 900]; drawnow set(f3, 'InvertHardcopy', 'off') if do_plotBoundaryOverlay f4 = figure(4); hold on f4t = tiledlayout('flow','TileSpacing','compact','Padding','compact'); boxsize = boxt1; for k = 1:nframes tmpbdy=bdy{k}; movieFrame=listFrames(k); %movieframe, usually same as im0:imL; % TIFF indices if ~do_fixedbox % adjust overall min, max after each frame minCol = min(minCol,min(tmpbdy(:,2)-dxadj(1))); maxCol = max(maxCol,max(tmpbdy(:,2)-dxadj(1))); minRow = min(minRow,min(tmpbdy(:,1)-dxadj(2))); maxRow = max(maxRow,max(tmpbdy(:,1)-dxadj(2))); boxsize = [minCol maxCol minRow maxRow]; end cm4=colormap(parula(fdelta)); % use as many colors as there frames thiscm4 = cm4(movieFrame-im0+1,:); nexttile(1) ax4=plot(tmpbdy(:,2)-dxadj(1),tmpbdy(:,1)-dxadj(2),'-','Color',thiscm4,'LineWidth',6); axis equal, hold on set(ax4.Parent,'Color',backcm) set(ax4.Parent,'Ydir','reverse') axis(boxsize) title(strcat('Boundary'),'FontSize',afs+12); ax4.Parent.CLim = [im0 imL+1]; %cc = colorbar('Ticks',im0+0.5:tickstep:imL+0.5,'TickLabels',{im0:tickstep:imL}); cc = colorbar('Ticks',im0-0.5+tickstep:tickstep:imL+0.5,'TickLabels',{im0-1+tickstep:tickstep:imL}); cc.Label.String = 'Frames'; end if ~do_fixedbox %extra padding when adjusting box for every GUV minCol = max( 1,minCol-5); maxCol = min(cols,maxCol+5); minRow = max( 1,minRow-5); maxRow = min(rows,maxRow+5); boxsize = [minCol maxCol minRow maxRow]; end hold off title(f4t,guvname,'FontSize',afs+14) ax4.Parent.FontName = 'Arial'; ax4.Parent.FontSize = afs+4; % ax4.Parent.XLim = boxsize(1:2); % ax4.Parent.YLim = boxsize(3:4); axis(boxsize) % width height f4.Position(3:4) = [500,425]; drawnow set(f4,'InvertHardcopy','off') % OUTPUTS figsout(1) = f3; figsout(2) = f4; if savefigs>0 print(f3,strcat(outputfolder,guvname,'_',tstamp1,'_fig_macrosummary'),'-dpng') print(f4,strcat(outputfolder,guvname,'_',tstamp1,'_fig_boundary'),'-dpng') end else % OUTPUTS figsout(1) = f3; if savefigs>0 print(f3,strcat(outputfolder,guvname,'_',tstamp1,'_fig_macrosummary'),'-dpng') end end end function figout = displayboundary(I1rgb,bdy,Borigs,Bkpoints) % this displays boundary f1= figure(1); clf(f1); tf1 = tiledlayout(2,4); tf1.Padding = 'compact'; tf1.TileSpacing = 'compact'; nexttile(1,[2 2]) imshow(I1rgb); hold on; plot(bdy(1,2),bdy(1,1),'go','LineWidth',10) plot(bdy(:,2),bdy(:,1),'g*','LineWidth',2) plot(Borigs(1,2),Borigs(1,1),'bo','LineWidth',10) plot(Borigs(:,2),Borigs(:,1),'b*','LineWidth',2) title('Raw boundary (blue)') hold off nexttile(3,[2 2]) imshow(I1rgb); hold on; plot(bdy(1,2),bdy(1,1),'og','LineWidth',10) plot(bdy(:,2),bdy(:,1),'*g','LineWidth',2) plot(Bkpoints(1,2),Bkpoints(1,1),'ob','LineWidth',10) plot(Bkpoints(:,2),Bkpoints(:,1),'*b','LineWidth',2) title('Spaced and smoothed boundary(blue)') hold off sgtitle('Adjustments to raw boundary (green) for circumcenter-based curvature computation') drawnow figout = f1; end function figout = displaymultipanel(fignum,fvars,I1rgb,I1pan,varargin) % recall: f2vars = v2struct(savgolay_smooth,computecurv,movieFrame,panelnames,panelcaxes,tmpbdy); v2struct(fvars); dmark = 45; FontSize = minfont; %min(minfont,floor(min(cols,rows)/10)); interpBdy = tmpbdy; npanel = 0.5*(nargin-2); f2= figure(fignum); clf(f2); tf2 = tiledlayout(2,2*npanel); tf2.Padding = 'compact'; tf2.TileSpacing = 'compact'; f2.Position(3:4) = [1500, 500]; Irgb = cell(npanel,1); Ivars = cell(npanel,1); Icax = cell(npanel,1); Irgb{1} = I1rgb; % background image Ivars{1} = I1pan.vars; % boundary quantities to overlay Icax{1} = I1pan.caxes; % color axis for boundary quantity if ~isempty(varargin) for i = 2:npanel Irgb{i} = varargin{2*(i-1)-1}; Ivars{i} = varargin{2*(i-1)}.caxes; Icax{i} = varargin{2*(i-1)}.caxes; end end for j = 1:npanel nexttile(2*j-1,[2 2]) imshow(Irgb{j}); hold on; f2ax1 = scatter3(interpBdy(:,2),interpBdy(:,1),zeros(length(interpBdy),1)... ,[],interpCurv,'filled','SizeData',dmark); if computecurv==1 title('Curvature (circum) along interpolated boundary') else title('Curvature (fitcircle) along interpolated boundary') end colorbar; colormap(parula); f2ax1.Parent.CLim = curvcax; %[-0.03,0.03]; %caxis([-0.11 0.11]); text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) hold off nexttile(3,[2 2]); % 5 if using flow imshow(I2rgb); hold on; f2ax2 = scatter3(interpBdy(:,2),interpBdy(:,1),zeros(length(interpBdy),1)... ,[],interpFluor,'filled','SizeData',dmark); title([fluorname(9:end),' normalized intensity along interpolated boundary']) colorbar; colormap(parula); f2ax2.Parent.CLim = tmpfluorcax; %[0,30];%caxis([0 41]); text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) hold off nexttile(5,[2 2]); % 9 if using flow imshow(I1rgb); hold on; if nargin == 7 interpVel = varargin{4}; f2ax3 = scatter3(interpBdy(:,2),interpBdy(:,1),zeros(length(interpBdy),1)... ,[],interpVel,'filled','SizeData',dmark); end title('Velocities along interpolated boundary') colorbar; colormap(parula); f2ax3.Parent.CLim = velcax; % [-5,5]; % caxis([-8.5 8.5]); text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) end hold off drawnow % adjust position f2.Position(3:4) = [1500, 500]; figout = f2; end % function figout = displayrest(I1rgb,I2rgb,f2vars,varargin) % v2struct(f2vars); % dmark = 45; % FontSize = minfont; %min(minfont,floor(min(cols,rows)/10)); % % if nargin >= 6 % f2= figure(2); clf(f2); % tf2 = tiledlayout(2,6); % tf2.Padding = 'compact'; % tf2.TileSpacing = 'compact'; % f2.Position(3:4) = [1500, 500]; % % interpBdy = varargin{1}; % interpCurv = varargin{2}; % interpGFP = varargin{3}; % % % % % nexttile(1,[2 2]) % imshow(I1rgb); % rows = size(I1rgb,1); cols = size(I1rgb,2); % hold on; % f2ax1 = scatter3(interpBdy(:,2),interpBdy(:,1),zeros(length(interpBdy),1)... % ,[],interpCurv,'filled','SizeData',dmark); % if computecurv==1 % title('Curvature (circum) along interpolated boundary') % else % title('Curvature (fitcircle) along interpolated boundary') % end % colorbar; colormap(parula); % f2ax1.Parent.CLim = curvcax; %[-0.03,0.03]; %caxis([-0.11 0.11]); % text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) % hold off % % nexttile(3,[2 2]); % 5 if using flow % imshow(I2rgb); % hold on; % f2ax2 = scatter3(interpBdy(:,2),interpBdy(:,1),zeros(length(interpBdy),1)... % ,[],interpGFP,'filled','SizeData',dmark); % title([fluorname(9:end),' normalized intensity along interpolated boundary']) % colorbar; colormap(parula); % f2ax2.Parent.CLim = tmpfluorcax; %[0,30];%caxis([0 41]); % text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) % hold off % % nexttile(5,[2 2]); % 9 if using flow % imshow(I1rgb); % hold on; % if nargin == 7 % interpVel = varargin{4}; % f2ax3 = scatter3(interpBdy(:,2),interpBdy(:,1),zeros(length(interpBdy),1)... % ,[],interpVel,'filled','SizeData',dmark); % end % title('Velocities along interpolated boundary') % colorbar; colormap(parula); % f2ax3.Parent.CLim = velcax; % [-5,5]; % caxis([-8.5 8.5]); % text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) % hold off % % drawnow % % adjust position % f2.Position(3:4) = [1500, 500]; % % figout = f2; % % else % f2= figure(2); clf(f2); % tf2 = tiledlayout(2,2); % tf2.Padding = 'compact'; % tf2.TileSpacing = 'compact'; % % % interpBdy = varargin{1}; % interpGFP = varargin{2}; % % nexttile(1,[2 2]) % imshow(I2rgb); % hold on; % f2ax2 = scatter3(interpBdy(:,2),interpBdy(:,1),zeros(length(interpBdy),1)... % ,[],interpGFP,'filled','SizeData',dmark); % title([fluorname(9:end),' normalized intensity along interpolated boundary']) % colorbar; colormap(parula); % f2ax2.Parent.CLim = tmpfluorcax; %[0,30];%caxis([0 41]); % text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) % hold off % % drawnow % % adjust position % f2.Position(3:4) = [500, 500]; % % figout = f2; % end % % end function varargout = interpboundary(bdy,n,varargin) % This function takes a boundary and, possibly, data and a number n and % interpolates the boundary and any data so that it is now of length n % for circular data: n = (# of desired unique points) + 1 % output is circular if input is circular itype = 'makima'; % avoids excessive undulations %itype = 'pchip'; if nargin~=nargout+1 disp('Error: Number of input and output data sets must equal') else d = diff(bdy); % difference between rows of bdy dist = sqrt(d(:,1).^2+d(:,2).^2); cumdist = [0;cumsum(dist)]; % (x) cumulative distances around boundary perim = max(cumdist); s = linspace(0,1,n)*perim; % (xq) n query points from 0 to perim % vq = interp1(x,v,xq,method) % interpolate (row,col) ie (y,x) coordinates across cumulative sum of boundary arc length row = interp1(cumdist,bdy(:,1),s,itype)'; col = interp1(cumdist,bdy(:,2),s,itype)'; B = [row col]; varargout = cell(length(varargin)+1,1); varargout{1}=B; for i=1:length(varargin) % if nans just return nans if all(isnan(varargin{i})) varargout{i+1} = nan(length(row),1); % disp('Error: Found NaNs in list') else % interpolate fluorescent signal,velocities,other quantities % across cumulative sum of boundary arc length varargout{i+1} = interp1(cumdist,varargin{i},s,itype)'; end end end end function varargout = aligninitboundary(I1rgb,bdy,z0,thetatop,varargin) % !!! make bdy quantities non-circular for alignment tmpbdy = bdy(1:end-1,:); %figure(71); imshow(I1rgb); hold on; axis equal; plot(tmpbdy(1,2),tmpbdy(1,1),'r*'); hold off tmpvarsin = cell(length(varargin),1); tmpvarsout = cell(length(varargin),1); for i=1:length(varargin) tmpvarsin{i} = varargin{i}(1:end-1); end x0 = z0(1); y0 = z0(2); theta = atan2(y0-tmpbdy(:,1), tmpbdy(:,2)-x0); % deltay = y0 - bdyrow; delta x = bdycol - x0; % theta here will be always between -pi to pi thetatop; % this is the angle pi away from thetacenter % thetatop will always be between -pi to pi %[~,idx]=(min(abs(theta-theta0))); [~,idx]=(min(abs(theta-thetatop))); % find index closest to thetapi % this will give index of point to align top of kymograph % m = 1-idx = negative number --> CW rotation of position 1 label; % shifts first m entries to end % theta0 marks top row of kymograph tmpbdyalign = circshift(tmpbdy, 1-idx); %figure(72); imshow(I1rgb); hold on; axis equal; plot(tmpbdyalign(1,2),tmpbdyalign(1,1),'r*'); hold off % make circular at end varargout{1} = [tmpbdyalign; tmpbdyalign(1,:)]; for i=1:length(varargin) tmpvarsout{i} = circshift(tmpvarsin{i},1-idx); % make circular at end varargout{i+1} = [tmpvarsout{i}; tmpvarsout{i}(1) ]; end end function varargout = alignboundary(newbdy,oldbdy,varargin) % newbdy and oldbdy are coordinates (y,x) relative to centroid % The analysis assumes that the new boundary is longer; if not, switch them % [alignedbdy,mindst,aligneddata1,...] = alignboundary(newbdy,oldbdy,data1,...) % !!! make bdy quantities non-circular for alignment tmpvarsin = cell(length(varargin),1); tmpvarsout = cell(length(varargin),1); for i=1:length(varargin) tmpvarsin{i} = varargin{i}(1:end-1); end sw = logical(length(newbdy)>length(oldbdy)); % make b1, b2 non-circular if sw %new bdy is longer- set this to b1 b1 = newbdy(1:end-1,:); b2 = oldbdy(1:end-1,:); else % old bdy is longer- set this to b1 b1 = oldbdy(1:end-1,:); b2 = newbdy(1:end-1,:); end n1 = length(b1); % number of points in longer bdy n2 = length(b2); % number of points in shorter bdy d=zeros(n2,1); % based on shorter bdy (b2) idx=[]; % add unique indices (between 1 and n1) to idx % until we get enough = n1-n2 while length(idx)<n1-n2 % difference in number of points deltan = n1-n2-length(idx); nidx=ceil(n1*rand(deltan,1)); % nidx = deltan random numbers between 1 and n1 idx=unique(sort([idx;nidx])); % idx = collect unique numbers in idx,nidx in ascending order end % find indices in 1:n1 but not in idx and sort % kidx has n1-(n1-n2) = n2 # of points % kidx collects n2 unique random indices between 1:n1, kidx = sort(setdiff(1:n1,idx))'; % b1s: list of n2 unique random boundary points from b1 % b1s has same length as b2 (shorter bdy) b1s = b1(kidx,:); dst = zeros(length(b1s),1); for i=0:length(b1s)-1 % find distance b/w x,y coords of b2 and b1s (both same # points) % cycle through all possible rotations % b1,b2 are listed CW initially d=b2-circshift(b1s,i);% rotate initial point in b1s CCW i positions dst(i+1)=mean(sqrt(sum(d'.^2))); end %figure(72) %plot(dst) [mindst,minidx]=min(dst); %minidx - 1 = how much b1s initial point should be rotated CCW (+circshift) if sw % new bdy longer, then newbdy = b1, so do same CCW direction as b1s newbdyalign = circshift(b1,minidx-1); varargout{1} = [newbdyalign; newbdyalign(1,:)]; % make circular for i=1:length(varargin) tmpvarsout{i} = circshift(tmpvarsin{i},minidx-1); varargout{i+1} = [tmpvarsout{i}; tmpvarsout{i}(1)]; end else % new bdy shorter, then newbdy = b2, so do opp CW but same amount as b1s % CCW direction for b1s = (minidx-1) % CW direction for b2 = - (minidx-1) newbdyalign = circshift(b2,-(minidx-1)); varargout{1} = [newbdyalign; newbdyalign(1,:)]; % make circular for i=1:length(varargin) tmpvarsout{i} = circshift(tmpvarsin{i},-(minidx-1)); varargout{i+1} = [tmpvarsout{i}; tmpvarsout{i}(1)]; end end end function varargout = alignboundary_ba(newbdy,oldbdy,varargin) % inputs: % newbdy, oldbdy, varargin: circular lists (repetition of first point at end of list) % varargin: vel, curv, fluor only % function determines alignment % create equal length boundary vectors (without repetition of first point) % rotate newbdy and compute sum of displacement map between the two ordered lists (newbdy and oldbdy) % aligned newbdy = rotation of newbdy that minimizes the sum of displacements % computes a motion tracking, with 1:1 correspondence % outputs: all aligned, circular, and interpolated bac to original newbdy length if newbdy was adjusted % varargout{1} = newbdy % varargout{2} = curvature for newbdy % varargout{3} = velocity (px/frame) computed based on pdist2 (NOT motion tracking) % varargout{4} = fluor for newbdy tmpvarsin = cell(1,length(varargin)); tmpvarsinu = cell(1,length(varargin)); tmpvarsout = cell(1,length(varargin)); tmpvarsoutC = cell(1,length(varargin)); varargout = cell(1,length(varargin)+1); % newbdy and oldbdy are circular lists = # unique points + 1 nnew = length(newbdy)-1; nold = length(oldbdy)-1; if nnew > nold % make oldbdy longer nlong = nnew; b1 = interpboundary(oldbdy,nlong+1); b2 = newbdy; % store additional arguments (keep circular) for i = 1:length(varargin) tmpvarsin{i} = varargin{i}; end elseif nnew < nold % make newbdy longer, and all other bdy quantities longer using interpolation nlong = nold; b1 = oldbdy; b2 = interpboundary(newbdy,nlong+1); for i = 1:length(varargin) [~,tmpvarsin{i}] = interpboundary(newbdy,nlong+1,varargin{i}); end else nlong = nnew; % same as nold bc lengths are equal in this case b1 = oldbdy; b2 = newbdy; % store additional arguments (keep circular) for i = 1:length(varargin) tmpvarsin{i} = varargin{i}; end end dst = zeros(nlong,1); % make bdy quantities non-circular (unique points only for interpolation and alignment) b1u = b1(1:end-1,:); b2u = b2(1:end-1,:); for i = 1:length(varargin) tmpvarsinu{i} = tmpvarsin{i}(1:end-1); end for i=0:nlong-1 % find distance b/w x,y coords of b2u and b1u (both same # points) % cycle through all possible rotations % b1u,b2u are listed CW initially d=b1u-circshift(b2u,i);% rotate initial point in b2u CCW i positions dst(i+1)=mean(sqrt(sum(d'.^2))); end % figure(72) % plot(dst) [~,minidx]=min(dst); %minidx - 1 = how much b2 initial point should be rotated CCW (+circshift) if nnew < nold % interpolation needed to downsample b2 back to original bdynew length b2align = circshift(b2u,minidx-1); % align b2 for i=1:length(varargin) tmpvarsout{i} = circshift(tmpvarsinu{i},minidx-1); %align other quantities end % make all quantities circular again before interpolation b2alignC = [b2align; b2align(1,:)]; for i=1:length(varargin) tmpvarsoutC{i} = [tmpvarsout{i}; tmpvarsout{i}(1)]; end % shorten b2alignC and other quants back to original newbdy length, circular output tmpvarsoutC_orig = cell(size(tmpvarsoutC)); [b2alignC_orig, tmpvarsoutC_orig{:}] = interpboundary(b2alignC,nnew+1,tmpvarsoutC{:}); varargout{1} = b2alignC_orig; for i = 1:length(varargin) varargout{1+i} = tmpvarsoutC_orig{i}; end else % nnew = nold OR nnew > nold % no interpolation needed bc b2align is already same length as bdynew b2align = circshift(b2u,minidx-1); varargout{1} = [b2align; b2align(1,:)]; % make circular for i=1:length(varargin) tmpvarsout{i} = circshift(tmpvarsinu{i},minidx-1); varargout{i+1} = [tmpvarsout{i}; tmpvarsout{i}(1)]; end end end function newvelT = velocitytrack(bdy1,c1,bdy0,c0,dt,I1rgb,pveltrack) % inputs: % newbdy, oldbdy: [row col], circular lists (repetition of first point at end of list) % newc, oldc: centroid positions [y0, x0] % pveltrack: structure array % computes velocity based on aligned boundaries (motion tracking) % create aligned/CENTERED equal length boundary vectors (without repetition of first point) % rotate newbdy and compute sum of displacement map between the two ordered lists (newbdy and oldbdy) % aligned newbdy = rotation of newbdy that minimizes the sum of displacements % computes net displacement = aligned disp + centroid disp % vel = net disp / dt % output: all aligned, circular, and interpolated back to original newbdy length if newbdy was adjusted % varargout{2} = velocity based on motion tracking v2struct(pveltrack); newbdy = bdy1 - c1; oldbdy = bdy0 - c0; % newbdy and oldbdy are circular lists = # unique points + 1 nnew = length(newbdy)-1; nold = length(oldbdy)-1; if nnew > nold % make oldbdy longer nlong = nnew; b1 = interpboundary(oldbdy,nlong+1); b2 = newbdy; elseif nnew < nold % make newbdy longer, and all other bdy quantities longer using interpolation nlong = nold; b1 = oldbdy; b2 = interpboundary(newbdy,nlong+1); else nlong = nnew; % same as nold bc lengths are equal in this case b1 = oldbdy; b2 = newbdy; end dst = zeros(nlong,1); % make bdy quantities non-circular (unique points only for interpolation and alignment) b1u = b1(1:end-1,:); b2u = b2(1:end-1,:); for i=0:nlong-1 % find distance b/w x,y coords of b2u and b1u (both same # points) % cycle through all possible rotations % b1u,b2u are listed CW initially d=b1u-circshift(b2u,i);% rotate initial point in b2u CCW i positions dst(i+1)=mean(sqrt(sum(d'.^2))); end % figure(72) % plot(dst) [~,minidx]=min(dst); %minidx - 1 = how much b2 initial point should be rotated CCW (+circshift) if nnew < nold % interpolation needed to downsample b2 back to original bdynew length b2ualign = circshift(b2u,minidx-1); % align b2 b2uorig = b2ualign + c1; % original coords (to account for centroid disp) b1uorig = b1u + c0; d0 = (b2uorig - b1uorig) ; % NET displacements between two non-circular lists (new - old bdy same length) disp = sqrt(sum(d0'.^2))'; %[in] = inpolygon(xq,yq,xv,yv) returns [in] indicating if the query points specified by xq and yq are inside or on the edge of the polygon area defined by xv and yv. % 1 if xq, yq is INSIDE or ON the polygon defined by the query points [in3] = inpolygon(b2uorig(:,2),b2uorig(:,1),b1uorig(:,2),b1uorig(:,1)); inout3 = 1-2*in3; svel3 = (1/dt)*disp.*inout3; % make all quantities circular again before interpolation b2uorigC = [b2uorig; b2uorig(1,:)]; b1uorigC = [b1uorig; b1uorig(1,:)]; velC = [svel3; svel3(1)]; % shorten b2alignC and dispC back to original newbdy length, circular output [b2_norig,velC_norig] = interpboundary(b2uorigC,nnew+1,velC); newvelT = velC_norig; else % nnew = nold OR nnew > nold % no interpolation needed bc b2align is already same length as bdynew b2ualign = circshift(b2u,minidx-1); b2uorig = b2ualign + c1; % original coords (to account for centroid disp) b1uorig = b1u + c0; d0 = (b2uorig - b1uorig) ; % NET displacements between two non-circular lists (new - old bdy same length) disp = sqrt(sum(d0'.^2))'; %[in] = inpolygon(xq,yq,xv,yv) returns [in] indicating if the query points specified by xq and yq are inside or on the edge of the polygon area defined by xv and yv. % 1 if xq, yq is INSIDE or ON the polygon defined by the query points [in3] = inpolygon(b2uorig(:,2),b2uorig(:,1),b1uorig(:,2),b1uorig(:,1)); inout3 = 1-2*in3; svel3 = (1/dt)*disp.*inout3; % make all quantities circular again and NO interpolation b2uorigC = [b2uorig; b2uorig(1,:)]; b1uorigC = [b1uorig; b1uorig(1,:)]; b2_norig = b2uorigC; velC = [svel3; svel3(1)]; newvelT = velC; end if do_plotVelseries cols = size(I1rgb,2); rows = size(I1rgb,1); npts = length(b2uorigC); idx = 1:round(npts/50):npts; % use list of boundary points that are equal and before any last step interpolation bdy = b2uorigC; obdy = b1uorigC; obdym = fliplr(obdy)'; % xcoords row1, ycoords row2 bdym = fliplr(bdy)'; % xcoords row1, ycoords row2 x12 = [obdym(1,idx); bdym(1,idx)]; y12 = [obdym(2,idx);bdym(2,idx)]; velCidx = velC(idx); f5 = figure(5);clf(f5); f5.Position(3:4) = [600,600]; tf5 = tiledlayout(2,2); tf5.TileSpacing = 'compact'; tf5.Padding='compact'; nexttile(1,[2,2]); imshow(0.1*I1rgb); axis equal; hold on; plot(bdy(:,2),bdy(:,1),'Color',[1 1 1],'LineWidth',5); plot(obdy(:,2),obdy(:,1),'Color',0.5*[1 1 1],'LineWidth',2,'LineStyle','--'); for i = 1:length(idx) %figure(5); f5ax1 = surf([x12(:,i) x12(:,i)],... [y12(:,i) y12(:,i)], ... zeros(2,2),velCidx(i)*ones(2,2),'FaceColor','none','EdgeColor','flat',... 'LineWidth',3); end colorbar; colormap(parula); f5ax1.Parent.CLim = [-5 5]; %caxis([-9, 9]); title('Boundary velocity (px/frame)') FontSize = minfont; %min(minfont,floor(min(cols,rows)/10)); text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize); hold off; drawnow f5.Position(3:4) = [600,600]; if savefigs>0 % create tiff series %pause(0.05) F5 = getframe(f5); if movieFrame == im0 imwrite(F5.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_vel.tif')); else imwrite(F5.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_vel.tif'),... 'WriteMode','append'); end end end % end of do_plotVelseries end function [I,Iorig]=readframe(fullFileName,frame) if ischar(fullFileName) % can read tif file Iorig = imread(fullFileName,frame); else % can input image directly Iorig = fullFileName(:,:,frame); end if ndims(Iorig)>2 Iorig = Iorig(:,:,2); end I = imgaussfilt(Iorig,1); I = medfilt2(I,[1 1]); end function [svel3,idx] = getvelocities(bdy,obdy,dt,pVel,varargin) % BSA v2struct(pVel); [dist, idx] = pdist2(obdy,bdy,'euclidean','Smallest',1); % dist contains K = 1 smallest pairwise distance to a point in obdy for each observation in bdy. % For each observation in Y, pdist2(X,Y,'smallest'/'largest',K) finds the K smallest or largest distances % by computing and comparing the distance values to all the observations in X. %idx contains: (row) index of point in X that was mapped to each element of Y. dist = dist'; % this option requires figure being opened %oROI1 = drawfreehand('Position', fliplr(obdy),'Visible','off'); %inout1 = 1-2*inROI(oROI1,bdy(:,2),bdy(:,1)); %oROI2 = images.roi.Freehand('Position',fliplr(obdy)); %inout2 = 1-2*inROI(oROI2,bdy(:,2),bdy(:,1)); %inROI: points cannot lie on the roi, only inside/outside. (see 'classify pixels that partially enclosed') % hold off % drawnow % vel = 1/dt * displacement (+ for outward, - for inward) % svel1 = (1/dt)*dist.*inout1; %[in] = inpolygon(xq,yq,xv,yv) returns [in] indicating if the query points specified by xq and yq are inside or on the edge of the polygon area defined by xv and yv. % 1 if xq, yq is INSIDE or ON the polygon defined by the query points [in3] = inpolygon(bdy(:,2),bdy(:,1),obdy(:,2),obdy(:,1)); inout3 = 1-2*in3; svel3 = (1/dt)*dist.*inout3; % figure(1) % imshow(I1rgb); hold on; axis equal; % plot(bdy(:,2),bdy(:,1),'r--'); % plot(obdy(:,2),obdy(:,1),'g--'); % plot(bdy(100,2),bdy(100,1),'ro'); % plot(obdy(100,2),obdy(100,1),'go'); % hold off; % drawnow %% plot displacements if do_plotVelseries I1rgb = varargin{1}; cols = size(I1rgb,2); rows = size(I1rgb,1); f5 = figure(5);clf(f5); obdym = fliplr(obdy(idx,:))'; % xcoords row1, ycoords row2 bdym = fliplr(bdy)'; % xcoords row1, ycoords row2 x12 = [obdym(1,:); bdym(1,:)]; y12 = [obdym(2,:);bdym(2,:)]; imshow(0.1*I1rgb); hold on; axis equal; plot(bdy(:,2),bdy(:,1),'Color',[1 1 1],'LineWidth',5); plot(obdy(:,2),obdy(:,1),'Color',0.5*[1 1 1],'LineWidth',2,'LineStyle','--'); for i = 1:length(bdy) % figure(5); f5ax = surf([x12(:,i) x12(:,i)],... [y12(:,i) y12(:,i)], ... zeros(2,2),svel3(i)*ones(2,2),'FaceColor','none','EdgeColor','flat',... 'LineWidth',3); end colorbar; colormap(parula); f5ax.Parent.CLim = [-5 5]; %caxis([-9, 9]); title('Boundary velocity (px/frame)') FontSize = min(minfont,floor(min(cols,rows)/10)); text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize); hold off; drawnow f5.Position(3:4) = [600,600]; if savefigs > 0 % create tiff series %pause(0.05) F5 = getframe(f5); if movieFrame == im0 imwrite(F5.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_vel.tif')); else imwrite(F5.cdata,strcat(outputfolder,guvname,'_',tstamp1,'_tseries_vel.tif'),... 'WriteMode','append'); end end end % end of do_plotVelseries end function [getFluor,getiBW,getoBW,getinner,getouter] = getIntensity(bdy,Iinit,BW,pfluor) % I will be a cell array based on fluorescent channels v2struct(pfluor); % extract parameters % if bits == 16 % maxfl = 65535; % end getFluor = cell(1,nfluor); %intensities getiBW = cell(1,nfluor); getoBW = cell(1,nfluor); getinner = cell(1,nfluor); getouter = cell(1,nfluor); count = 0; for k = 1:nchannels if ismember(k,fluorch) count = count + 1; del = 4;% 4 originally # of points left or right used to find perpendicular line segment lperp = nsize*8;% px outer and px inner for perpendicular line segment pin = nsize*8; % px to erode for luminal mask pout1 = nsize*8; % px to dilate for 1st outer ring of external mask pout2 = nsize*16; % px to dilate for 2nd outer ring of external mask ptop = [1]; % prop of top gfp intensities to use along lperp % KEEP ptop = 1; USE PLPERP to adjust how far outward to compute top gfp intensities npts = length(bdy); fluor = zeros(npts,1); % for ptop1 linner = zeros(npts,2); louter = zeros(npts,2); iBW = imerode(BW,strel('disk',round(pin),4));% 6,4 oBW1 = imdilate(BW,strel('disk',round(pout1),4));% 6,4 oBW2 = imdilate(BW,strel('disk',round(pout2),4));% 16,4 oBW = immultiply(imcomplement(oBW1),oBW2); Ifl = Iinit{k}; % single channel fluorescent signal Iflrgb = repmat(imadjust(Ifl),1,1,3); cols = size(Ifl,2); rows = size(Ifl,1); cent = centroid; for pt=1:npts x0 = bdy(pt,2); y0 = bdy(pt,1); cb = circshift(bdy,del+1-pt); seg = cb(1:2*del+1,[2 1]); [z, r, ~] = fitcircle(seg); % z center, r radius if(isnan(r)) r=rold; z=zold; end xc = z(1); yc = z(2); % find normal to curve and extend line segment through in polar coords theta = atan2(yc-y0,x0-xc); % dy = yc - row# ; dx = col# - xc d = -lperp:0.5:lperp;% 0.5 increment originally x = x0+d*cos(theta); % x(1) = x0 - lperp*cosA ---> where x0 is the boundary point y = y0-d*sin(theta);% because row coord increases going down [~,idx]= min( [norm([x(1),y(1)]-cent), norm([x(end),y(end)]-cent)] ); lend = round(plperp*length(x)); % end point of line segment to consider if idx == 1 xord = x(1:lend); yord= y(1:lend); linner(pt,:)=[xord(1),yord(1)]; % first row of inner is coords of innermost point of line segment (d=-lperp) louter(pt,:)=[xord(end),yord(end)]; % first row of outer is coords of outermost point of line segment (d=+lperp) else % idx == 2 xa = fliplr(x); ya = fliplr(y); xord = xa(1:lend)'; yord = ya(1:lend)'; linner(pt,:)=[xord(1),yord(1)]; % first row of inner is coords of innermost point of line segment (d=-lperp) louter(pt,:)=[xord(end),yord(end)]; % first row of outer is coords of outermost point of line segment (d=+lperp) end rold = r; zold = z; tmptmp = improfile(Ifl,xord,yord,'bilinear'); %figure(33); imshow(imadjust(I)); hold on; plot(x,y,'r','LineWidth',5); hold off; %drawnow tmp = sort(tmptmp,'descend','MissingPlacement','last'); if do_top3px % mean of just top 3 pixels (less sensitive to spatial thickening) fluor(pt) = mean(tmp(1:3)); else % mean fluor(pt) = mean(tmp(1:ptop*end)); % top percentage to use % median %fluor(pt) = double(median(tmp(1:ptop*end))); % top percentage to use end end % membrane gfp normalized to background and internal signal % (gfp_memb - mean(background)) / (mean(gfp_internal) - mean(background)) %gfp = (gfp-mean(I(oBW)))/(mean(I(iBW))-mean(I(oBW))); % sort and take mean of bottom fraction of cytosolic (ibw) and external (obw) tmpibw = sort(Ifl(iBW),'MissingPlacement','last');% ascending list of gfp values in mask iBW tmpobw = sort(Ifl(oBW),'MissingPlacement','last');% ascending list of gfp values in mask oBW mibw = mean(tmpibw(1:round(0.2*end))); mobw = mean(tmpobw(1:round(0.2*end))); %mibw = double(median(tmpibw(1:round(0.2*end)))); %mobw = double(median(tmpobw(1:round(0.2*end)))); % if k > 1 % fluor = (fluor-mobw)/(mibw-mobw); % elseif k == 1 % fluor = (fluor-mobw)/maxfl; % end %% NO LUMEN division here: % fluor = (fluor-mobw)/(mibw-mobw); fluor = (fluor-mobw); %% getoBW{count} = oBW; getiBW{count} = iBW; getFluor{count} = fluor; getinner{count} = linner; getouter{count} = louter; if k > 1 && do_membnorm getFluor{count} = fluor./getFluor{1}; % divide by membrane channel end end % if ismember k fluorch end %end of nfluors end function [allkymos, alignall,varargout] = makekymos(dataRaw,kf) % allkymos : cell array, 1 by nbdyquants % each entry is a matrix of the kymograph to be plotted % rows = nbins or nframes, cols = 6 (vel, curv, dcumul, 3 fluor channels % alignall : cell array, nframes by (nbdyquants + 1) % each entry is a vector of aligned boundary quantities ... % (except entries in column 1: boundary points have 2 columns) % rows = nboundarypts or nbins or nframes, % cols = 1 usually, except boundary points which have 2 % alignall_long : same as alignall, but using full nboundarypoints always % this will be varargout, added in case where binning occurs % kymofixedpts = dataRaw.kymofixedpts; % kymofixedtheta = dataRaw.kymofixedtheta; guvname = dataRaw.guvname; bdyall = dataRaw.allBoundaries; % cell array nframes by 1 boundary points bdycurv = dataRaw.allCurvatures; % cell array nframes by 1 bdyvel = dataRaw.allpdist2Velocities; % cell array nframes by 1 bdypd2dcumul = dataRaw.allpdist2Dcumuls; bdyperim = dataRaw.allPerims; fluorall = dataRaw.allIntensities; % cell array nframes by nfluor cent = dataRaw.allCentroids; % matrix: nframes by 2 leng = dataRaw.allLengths; % vector: nframes by 1 im0 = dataRaw.im0; imL = dataRaw.imL; nframes = dataRaw.nframes; % note nframes does not always equal imL+1-im0 nfluor = dataRaw.nfluor; % number of fluorescent channels membch = dataRaw.membch; fluorch = dataRaw.fluorch; listFrames = dataRaw.listFrames; % may not be consecutive frames during testing thetacenter = dataRaw.thetacenter; % angle for center of kymograph fnames = dataRaw.fnames; savefigs = dataRaw.savefigs; outputfolder = dataRaw.outputfolder; tstamp1 = dataRaw.tstamp1; do_plotVelseries = dataRaw.do_plotVelseries; nbins = dataRaw.nbins; pvary = dataRaw.pvary; allIrgb = dataRaw.allIrgb; veltrack1 = dataRaw.veltrack1; % fixedkymo: 0 for pdist2, 1 for motiontracking veltrack2 = dataRaw.veltrack2; % variable kymo: 0 for pdist2, 1 for motiontracking dt = dataRaw.dt; deltaf = dataRaw.deltaf; do_movavg = dataRaw.do_movavg; s_movavg = dataRaw.s_movavg; minfont = dataRaw.minfont; % if no smoothing average, set s to 0 if ~do_movavg s_movavg = 0; end % new variables alignbdy = cell(nframes,1); alignvel = cell(nframes,1); aligncurv = cell(nframes,1); alignfluor = cell(nframes,nfluor); aligndcumul = cell(nframes,1); nbdyquants = 3+nfluor; % vel, curv, dcumul kymos + fluor kymos %alignall_long = cell(nframes,nbdyquants+1); % col: bdypositions, vel, curv, dcumul, fluor kymos if kf > 0 %kymofixedpts || kymofixedtheta n = max(leng); % max # of boundary points in circular list over all frames alignangles = cell(nframes,1); alignDtheta = cell(nframes,1); alignbinidx = cell(nframes,1); % bdy point bin indices for every frame if kf ==2 %kymofixedtheta edges = - pi : 2*pi/nbins : pi ; else edges = 1 : (n-1)/nbins : n; end varargout{1} = cell(nframes,3); %for alignangles, alignDtheta, alignbinidx varargout{2} = edges'; else nvary = round(pvary*bdyperim); end %-------------------------- for j=1:nframes movieFrame = listFrames(j); % snake interpolation if having fixed # of points or theta bins in kymograph if kf > 0 %kymofixedpts || kymofixedtheta fluor = cell(1,nfluor); % pre-allocation is required % interpolate bdy quantities % pass fluor as comma separated list, not as a cell array [tmp,vel,curv,pd2dcumul,fluor{:}] = interpboundary(bdyall{j},n,bdyvel{j},... bdycurv{j},bdypd2dcumul{j},fluorall{j,:}); % interpBdy and interpCurv are circular; others just very close curv(end) = curv(1); vel(end) = vel(1); pd2dcumul(end) = pd2dcumul(1); % make list of fluor circular and keep as cell array matfluor = cell2mat(fluor); matfluor(end,:) = matfluor(1,:); fluor = mat2cell(matfluor,length(fluor{1}),ones(1,nfluor)); % current centroid x0 = cent(j,1); % centroids are (col, row), (x, y) y0 = cent(j,2); else fluor = cell(1,nfluor); % pre-allocation is required % interpolate bdy quantities [tmp,vel,curv,pd2dcumul,fluor{:}] = interpboundary(bdyall{j},nvary(j),bdyvel{j},... bdycurv{j},bdypd2dcumul{j},fluorall{j,:}); % interpBdy and interpCurv are circular; others just very close curv(end) = curv(1); vel(end) = vel(1); pd2dcumul(end) = pd2dcumul(1); % make list of fluor circular and keep as cell array matfluor = cell2mat(fluor); matfluor(end,:) = matfluor(1,:); fluor = mat2cell(matfluor,length(fluor{1}),ones(1,nfluor)); % current centroid x0 = cent(j,1); % centroids are (col, row), (x, y) y0 = cent(j,2); % tmp = bdyall{j}; % (row col) (y, x) % x0 = cent(j,1); % centroids are (col, row), (x, y) % y0 = cent(j,2); % vel = bdyvel{j}; % curv = bdycurv{j}; % pd2dcumul = bdypd2dcumul{j}; % fluor = fluorall(j,:); %cell array end %-------------------------------------------------------------------- % We now align the boundaries. If this is the first frame, we realign % it so that the angle 'theta0' corresponds roughly to the point % that is at the center of the kymograph (assuming symmetric changes). % If this is not the first frame, we align it so that the points are closest % to the previous frame, relative to the centroid. %-------------------------------------------------------------------- % thetapi = angle (rad) of initial point in aligned list, will be pi away from thetacenter thetastar = wrapToPi(deg2rad(thetacenter)); % reference angle, center of kymograph thetapi= wrapToPi(thetastar+pi); % top of kymograph, % make output from top to bottom y-axis: +pi to 0 to -pi (boundary points always go clockwise) newfluor = cell(size(fluor)); % pre-allocate is required here if j==1 [newtmp,newvel,newcurv,newpd2dcumul,newfluor{:}] = aligninitboundary(allIrgb{j,1},tmp,[x0 y0],thetapi,vel,curv,pd2dcumul,fluor{:}); if kf > 0 %kymofixedpts || kymofixedtheta newdcumul = zeros(size(newvel)); else % kymo varying, initial frame always 0 newdcumul = newpd2dcumul; end else % after first frame, align/track based on previous frame % newbdy: (y, x) relative to new centroid % oldbdy: (y, x) relative to old centroid [newtmp,newvel,newcurv,newpd2dcumul,newfluor{:}] = alignboundary_ba(tmp-[y0 x0],oldtmp-[y0old x0old],vel,curv,pd2dcumul,fluor{:}); % if veltrack % pveltrack=v2struct(movieFrame,do_plotVelseries,im0,savefigs,outputfolder,tstamp1,guvname,minfont); % [newvel] = velocitytrack(tmp,[y0 x0],oldtmp,[y0old x0old],dt,allIrgb{j,membch},pveltrack); % end if kf > 0 && veltrack1 %kymofixedpts || kymofixedtheta pveltrack=v2struct(movieFrame,do_plotVelseries,im0,savefigs,outputfolder,tstamp1,guvname,minfont); [newvel] = velocitytrack(tmp,[y0 x0],oldtmp,[y0old x0old],dt,allIrgb{j,membch},pveltrack); newdcumul = olddcumul + dt*newvel; elseif kf > 0 && ~veltrack1 newdcumul = olddcumul + dt*newvel; elseif kf == 0 && veltrack2 % kymo varying use pmotion tracking for vel pveltrack=v2struct(movieFrame,do_plotVelseries,im0,savefigs,outputfolder,tstamp1,guvname,minfont); [newvel] = velocitytrack(tmp,[y0 x0],oldtmp,[y0old x0old],dt,allIrgb{j,membch},pveltrack); % can only use pdist2 for dcumul newdcumul = newpd2dcumul; elseif kf == 0 && ~veltrack2 % kymovarying and no motiontrack for vel % use pdist2 for both velocity and dcumul % newvel already computed above newdcumul = newpd2dcumul; end newtmp = newtmp+[y0 x0]; end x0old = x0; y0old = y0; if kf > 0 % kymofixedpts || kymofixedtheta olddcumul = newdcumul; end %-------------------------------------------------------------------- % We are ready to save this frame's boundary. - angle computations, moving average %-------------------------------------------------------------------- % these are all circular lists if kf > 0 %kymofixedpts || kymofixedtheta % compute angles for every boundary point in frame j % compute displacement relative to thetastar (center of kymograph) if kf == 2 %kymofixedtheta %figure(10); imshow(allIrgb{1,1},[]); hold on; plot(tmp(1:50:200,2),tmp(1:50:200,1),'ro');plot(tmp(1,2),tmp(1,1),'g*');hold off; drawnow %figure(11); imshow(allIrgb{1,1},[]); hold on; plot(newtmp(1:50:200,2),newtmp(1:50:200,1),'ro');plot(newtmp(1,2),newtmp(1,1),'g*');hold off; drawnow % deltay = y0 - bdyrow; delta x = bdycol - x0; angles = atan2(y0-newtmp(:,1),newtmp(:,2)-x0); % tmp [row, col] D = wrapToPi(angles - thetastar); % angles relative to thetastar alignangles{j} = angles; % raw angles from -pi to pi alignDtheta{j} = D; % angles relative to thetastar else D = 1:n; %indices of points will be binned based on edges defined above for kymofixedpts end % if no moving average, s = 0; snewfluor = cell(size(newfluor)); [snewvel, snewcurv,snewdcumul,snewfluor{:}] = mavgsmoothboundary(s_movavg,newvel,newcurv,newdcumul,newfluor{:}); %alignbdy{j} = newtmp; binfluor = cell(1,nfluor); %[binidx, binvel, bincurv, bindcumul,binfluor{:}] = binbdyquant(kymofixed,D,edges,snewvel,snewcurv,snewdcumul,snewfluor{:}); [binidx, bintmprow, bintmpcol, binvel, bincurv, bindcumul,binfluor{:}] = ... binbdyquant(kf,D,edges,newtmp(:,1),newtmp(:,2),snewvel,snewcurv,snewdcumul,snewfluor{:}); alignbdy{j} = [bintmprow,bintmpcol]; % ============ for kymofixedtheta or fixedpts % we want to align quantities to be centered around a constant thetastar - done using D % for kymofixedpts % just want to bin quantities in fewer points alignbinidx{j} = binidx; alignvel{j} = binvel; aligncurv{j} = bincurv; aligndcumul{j} = bindcumul; alignfluor(j,:) = binfluor; % cell array, row: frame j, col: fluor channels testfornans = any(isnan(bincurv)); if testfornans kf j disp('some NaNs are present'); %return end oldtmp = newtmp; else % # of points in kymograph columns vary as perimeter changes % if no moving average, s = 0; snewfluor = cell(size(newfluor)); [snewvel, snewcurv,snewdcumul,snewfluor{:}] = mavgsmoothboundary(s_movavg,newvel,newcurv,newdcumul,newfluor{:}); alignbdy{j} = newtmp; alignvel{j} = snewvel; aligncurv{j} = snewcurv; aligndcumul{j} = snewdcumul; alignfluor(j,:) = snewfluor; % cell array, row: frame j, col: fluor channels oldtmp = newtmp; %alignall_long(j,[1:3,5:nbdyquants]) = [{newtmp},{snewvel},{snewcurv},snewfluor]; end end % end of nframes %compute cumulative displacement if kf > 0 %kymofixedtheta || kymofixedpts % v1 = cell2mat(alignvel'); % d1 = dt*v1; % dcumul = zeros(size(v1)); % dcumul(:,2:end) = cumsum(d1(:,2:end),2); % aligndcumul = mat2cell(dcumul,[nbins],ones(1,nframes))'; % transpose to match column format varargout{1} = [alignangles,alignDtheta,alignbinidx]; end %alignall_long(:,4) = aligndcumul; % will be empty for when usual kymograph where npoints varies allkymos = cell(1,nbdyquants); % every element is a matrix (max #boundary points by nframes) alignall = cell(nframes,nbdyquants+1); % aligned vel, curv, dcumul + fluor % (every element is different size matrix as # boundary points change) alignall(:,1:4) = [alignbdy,alignvel, aligncurv, aligndcumul]; alignall(:,5:end) = alignfluor; % ignore repeat for circ lists and circshift being used later for smoothing ky = cell(1,nbdyquants); % temporary variable for frame j for i = 1:nbdyquants % use deltaf to include blank columns in kymograph for skipped frames if kf > 0 %kymofixedtheta || kymofixedpts allkymos{i} = nan(nbins,deltaf); lkym = nbins*ones(nframes,1); % kymograph rows = nbins , fixed ell = zeros(nframes,1); % no padding needed else % using max(leng) will plot the repeat at end of circular lists % keeping it circular makes no difference bc so hard to see repeat allkymos{i} = nan(max(nvary),deltaf); lkym = nvary; % kymograph rows varies frame to frame % ell: half of difference between current bdy length and max bdy length ell = floor((max(nvary)-nvary)/2); % adjust every frame end for j = 1:nframes % frame j, i+1 gives bdyquantity of interest: keep circ fluor list; or all bins if doing binning ky{i} = alignall{j,i+1}; % use ell to center data jj = listFrames(j)-im0+1; % column index accounting for skipped frames allkymos{i}((ell(j)+1):(ell(j)+lkym(j)),jj)=ky{i}+0; end end % end of nbdyquants end % function function varargout = binbdyquant(kf,Dtheta,edges,varargin) nbins = length(edges)-1; Y = discretize(Dtheta,edges); % assign bin indices for every Dtheta, fixed for a given frame nvar = length(varargin); varargout = cell(1,nvar+1); if kf == 2 %bins go from -pi to pi (based on edges); % shift to be indexed from +pi to -pi Y1 = nbins-Y+1; % assign bin indices in reverse order varargout{1} = Y1; else Y1 = Y; varargout{1} = Y1; end for i = 1:nvar binX = nan(nbins,1); X = varargin{i}; for j = 1:nbins % find all data points assigned to bin index j binidx = (Y1==j); % same for every varargin, but changes placement of 1s as bins change % assign mean of data points in binidx to proper row in binned data binX(j) = mean(X(binidx)); end varargout{i+1} = binX; end end function varargout = v2struct(varargin) %% Description % v2struct has dual functionality in packing & unpacking variables into structures and % vice versa, according to the syntax and inputs. % % Function features: % * Pack variables to structure with enhanced field naming % * Pack and update variables in existing structure % * Unpack variables from structure with enhanced variable naming % * Unpack only specific fields in a structure to variables % * Unpack without over writing existing variables in work space % % In addition to the obvious usage, this function could by highly useful for example in % working with a function with multiple inputs. Packing variables before the call to % the function, and unpacking it in the beginning of the function will make the function % call shorter, more readable, and you would not have to worry about arguments order any % more. Moreover you could leave the function as it is and you could pass same inputs to % multiple functions, each of which will use its designated arguments placed in the % structure. % %- Syntax % Pack % S = v2struct % S = v2struct(x,y,z,...) % S = v2struct(fieldNames) % S = v2struct(A,B,C,..., fieldNames) % S = v2struct(x,..., nameOfStruct2Update, fieldNames) % v2struct % v2struct(x,y,z,...) % v2struct(fieldNames) % v2struct(A,B,C,..., fieldNames) % v2struct(x,..., nameOfStruct2Update, fieldNames) % % Unpack % v2struct(S) % [a,b,c,...] = v2struct(S) % v2struct(S,fieldNames) % [a,b,c,...] = v2struct(S,fieldNames) % %- Inputs & Outputs % Pack - inputs % x,y,z,... - any variable to pack. can be replaced by fieldNames below. % nameOfStruct2Update - optional, name of structure to update if desired. % fieldNames - optional, cell array of strings, which must include a cell % with the string 'fieldNames' and must be the last input. % Pack - outputs % S - the packed structure. If there is no output argument then a structure named % Sv2struct would be created in the caller workspace. % % Unpack - inputs % S - name of structure to be unpacked. % fieldNames - optional, cell array of strings, which must include a cell with the % string 'fieldNames' and must be the last input. % Unpack - outputs % a,b,c,... - variables upacked from the structure. % if there are no output arguments then variables would be created in % the caller workspace with naming according to name of inputs. % %- Examples % % see 'Usage example' section below for convenient presentation of these examples. % % % NOTE: whenever using filedNames cell array please note the following % % 1. fieldNames cell array must include a cell with the string 'fieldNames' % % 2. fieldNames cell array input must be the last input. % % % Pack % x = zeros(3); x2 = ones(3); y = 'Testing123'; z = cell(2,3); % fieldNames1 = {'fieldNames', 'x', 'y', 'z'}; % fieldNames2 = {'fieldNames', 'a', 'b', 'c'}; % fieldNames3 = {'fieldNames', 'x'}; % nameOfStruct2Update = 'S'; % % % The four examples below return structure S with same values however the % % structure's field names are defined differently in every syntax. % % Example 1. % % structure field names defined by variables names. % S = v2struct(x,y,z) % % Example 2. % % structure field names defined according to the cell array fieldNames. % % NOTE: variables with the names in fieldNames1 must exist in the caller workspace. % S = v2struct(fieldNames1) % % Example 3. % % same as #1. but arguments are passed explicitly % S = v2struct(zeros(3), 'Testing123', cell(2,3), fieldNames1) % % Example 4. % % field names defined by content of fieldNames2 while % % the values are set according to the passed arguments. In this case the structure % % S returned would be: S.a=x, S.b=y, S.c=z % S = v2struct(x,y,z, fieldNames2) % % % Example 5. % % update structure S. The fields that would be updated are according to content % % of fieldNames3. Note that you must pass a variable with the name % % 'nameOfStruct2Update' placed before 'fieldNames3'. This variable should contain % % the name of the structure you want to update as a string. Also note that if you % % set an output structure name which is different than the one stated in % % nameOfStruct2Update a new structure would be created and the structure that was % % meant to be updated would not get updated. % S.oldField = 'field to be saved for future use' % S = v2struct(x2, nameOfStruct2Update, fieldNames3) % % % Example 6. % % pack all variables in caller workspace. Call without input arguments. % S = v2struct % % % The following examples return the same results as the examples above but the % % structure would be returned with the default name 'Sv2struct'. Be cautious as % % this might lead to overriding of arguments. % % Example 7. % v2struct(x,y,z) % % Example 8. % v2struct(fieldNames1) % % Example 9. % v2struct(zeros(3), 'Testing123', cell(2,3), fieldNames1) % % Example 10. % v2struct(x,y,z, fieldNames2) % % Example 11. % S.oldField = 'field to be saved for future use' % v2struct(x2, nameOfStruct2Update, fieldNames3) % % Example 12. % v2struct % % % Unpack % clear S x x2 y z fieldNames1 fieldNames2 fieldNames3 nameOfStruct2Update % S.x = zeros(3); S.y = 'Testing123'; S.z = cell(2,3); % fieldNames3 = {'fieldNames','x','z'}; % % % Example 1. % % This example creates or overwrites variables x, y, z in the caller with the % % contents of the corresponding named fields. % v2struct(S) % % % Example 2. % % This example assigns the contents of the fields of the scalar structure % % S to the variables a,b,c rather than overwriting variables in the caller. If % % there are fewer output variables than there are fields in S, the remaining fields % % are not extracted. % [a,b,c] = v2struct(S) % % % Example 3. % % This example creates or overwrites variables x and z in the caller with the % % contents of the corresponding named fields. % v2struct(S, fieldNames3) % % % Example 4. % % This example assigns the contents of the fields 'x' and 'z' defined by % % fieldNames3 of the scalar structure S to the variables a and b rather than % % overwriting variables in the caller. If there are fewer output variables than % % there are fields in S, the remaining fields are not extracted. % [a,b] = v2struct(S, fieldNames3) % % % This example unpacks variables 'y' and 'z' only without overwriting variable 'x'. % % NOTE the addition of the field named 'avoidOverWrite' to the structure to be % % unpacked. This is mandatory in order to make this functionality work. The % % contents of this field can be anything, it does not matter. % S.avoidOverWrite = ''; % x = 'do not overwrite me'; % v2struct(S) % %- Usage example (includes sub-functions) % 1. run attached v2structDemo1.m file for on screen presentation of examples. % 2. run attached v2structDemo2.m file and read comments in file for a suggestion of % how to use v2struct in managing input to other functions with improved usability. % %- Revision history % 2011-05-19, Adi N., Creation % 2011-05-29, Adi N., Update structure added, some documentation and demo function changes % 2011-06-02, Adi N., Fixed updating structure functionality % 2011-06-05, Adi N., Added functionality: avoid overwritring existing variables, added % unpacking examples to demo1 .m file. % 2011-06-30, Adi N., fieldNames usage corrected, now must include a specific string to % be triggered. Documentation enhanced. Code tweaked. % 2011-07-14, Adi N., Fixed bug in packing with variables only % 2011-08-14, Adi N., Clarified warning and error when packing/unpacking with % fieldNames. % 2011-09-12, Adi N., Added easy packing of all variables in caller workspace (thanks % to Vesa Lehtinen for the suggestion), fixed bug in warning % handling in packing case, edited comments. % % Inspired by the function: mmv2truct - D.C. Hanselman, University of Maine, Orono, ME % 04469 4/28/99, 9/29/99, renamed 10/19/99 Mastering MATLAB 5, Prentice Hall, % ISBN 0-13-858366-8 % parse input for field names if isempty(varargin) gotCellArrayOfStrings = false; toUnpackRegular = false; toUnpackFieldNames = false; gotFieldNames = false; else gotCellArrayOfStrings = iscellstr(varargin{end}); toUnpackRegular = (nargin == 1) && isstruct(varargin{1}); if toUnpackRegular fieldNames = fieldnames(varargin{1})'; nFields = length(fieldNames); end gotFieldNames = gotCellArrayOfStrings & any(strcmpi(varargin{end},'fieldNames')); if gotFieldNames fieldNamesRaw = varargin{end}; % indices of cells with actual field names, leaving out the index to 'fieldNames' cell. indFieldNames = ~strcmpi(fieldNamesRaw,'fieldNames'); fieldNames = fieldNamesRaw(indFieldNames); nFields = length(fieldNames); end toUnpackFieldNames = (nargin == 2) && isstruct(varargin{1}) && gotFieldNames; end % Unpack if toUnpackRegular || toUnpackFieldNames struct = varargin{1}; assert(isequal(length(struct),1) , 'Single input nust be a scalar structure.'); CallerWS = evalin('caller','whos'); % arguments in caller work space % update fieldNames according to 'avoidOverWrite' flag field. if isfield(struct,'avoidOverWrite') indFieldNames = ~ismember(fieldNames,{CallerWS(:).name,'avoidOverWrite'}); fieldNames = fieldNames(indFieldNames); nFields = length(fieldNames); end if toUnpackRegular % Unpack with regular fields order if nargout == 0 % assign in caller for iField = 1:nFields assignin('caller',fieldNames{iField},struct.(fieldNames{iField})); end else % dump into variables for iField = 1:nargout varargout{iField} = struct.(fieldNames{iField}); end end elseif toUnpackFieldNames % Unpack with fields according to fieldNames if nargout == 0 % assign in caller, by comparing fields to fieldNames for iField = 1:nFields assignin('caller',fieldNames{iField},struct.(fieldNames{iField})); end else % dump into variables assert( isequal(nFields, nargout) , ['Number of output arguments',... ' does not match number of field names in cell array']); for iField = 1:nFields varargout{iField} = struct.(fieldNames{iField}); end end end % Pack else % build cell array of input names CallerWS = evalin('caller','whos'); inputNames = cell(1,nargin); for iArgin = 1:nargin inputNames{iArgin} = inputname(iArgin); end nInputs = length(inputNames); % look for 'nameOfStruct2Update' variable and get the structure name if ~any(strcmpi(inputNames,'nameOfStruct2Update')) % no nameOfStruct2Update nameStructArgFound = false; validVarargin = varargin; else % nameOfStruct2Update found nameStructArgFound = true; nameStructArgLoc = strcmp(inputNames,'nameOfStruct2Update'); nameOfStruct2Update = varargin{nameStructArgLoc}; % valid varargin with just the inputs to pack and fieldNames if exists validVarargin = varargin(~strcmpi(inputNames,'nameOfStruct2Update')); % valid inputNames with just the inputs name to pack and fieldNames if exists inputNames = inputNames(~strcmpi(inputNames,'nameOfStruct2Update')); nInputs = length(inputNames); % copy structure from caller workspace to enable its updating if ismember(nameOfStruct2Update,{CallerWS(:).name}) % verify existance S = evalin('caller',nameOfStruct2Update); else error(['Bad input. Structure named ''',nameOfStruct2Update,... ''' was not found in workspace']) end end % when there is no input or the input is only variables and perhaps % also nameOfStruct2Update if ~gotFieldNames % no input, pack all of variables in caller workspace if isequal(nInputs, 0) for iVar = 1:length(CallerWS) S.(CallerWS(iVar).name) = evalin('caller',CallerWS(iVar).name); end % got input, check input names and pack else for iInput = 1:nInputs if gotCellArrayOfStrings % called with a cell array of strings errMsg = sprintf(['Bad input in cell array of strings.'... '\nIf you want to pack (or unpack) using this cell array as'... ' designated names'... '\nof the structure''s fields, add a cell with the string'... ' ''fieldNames'' to it.']); else errMsg = sprintf(['Bad input in argument no. ', int2str(iArgin),... ' - explicit argument.\n'... 'Explicit arguments can only be called along with a matching'... '\n''fieldNames'' cell array of strings.']); end assert( ~isempty(inputNames{iInput}), errMsg); S.(inputNames{iInput}) = validVarargin{iInput}; end % issue warning for possible wrong usage when packing with an input of cell array of % strings as the last input without it containing the string 'fieldNames'. if gotCellArrayOfStrings name = inputNames{end}; % input contains structure and a cell array of strings if (nargin == 2) && isstruct(varargin{1}) msgStr = [inputNames{1},''' and ''',inputNames{2},''' were']; % input contains any arguments with an implicit cell array of strings else msgStr = [name, ''' was']; end warnMsg = ['V2STRUCT - ''%s packed in the structure.'... '\nTo avoid this warning do not put ''%s'' as last v2struct input.'... '\nIf you want to pack (or unpack) using ''%s'' as designated names'... ' of the'... '\nstructure''s fields, add a cell with the string ''fieldNames'' to'... ' ''%s''.']; fprintf('\n') warning('MATLAB:V2STRUCT:cellArrayOfStringNotFieldNames',warnMsg,msgStr,... name,name,name) end end % fieldNames cell array exists in input elseif gotFieldNames nVarToPack = length(varargin)-1-double(nameStructArgFound); if nVarToPack == 0 % no variables to pack for iField = 1:nFields S.(fieldNames{iField}) = evalin('caller',fieldNames{iField}); end % else - variables to pack exist % check for correct number of fields vs. variables to pack elseif ~isequal(nFields,nVarToPack) error(['Bad input. Number of strings in fieldNames does not match',... 'number of input arguments for packing.']) else for iField = 1:nFields S.(fieldNames{iField}) = validVarargin{iField}; end end end % if ~gotFieldNames if nargout == 0 assignin( 'caller', 'Sv2struct',S ); else varargout{1} = S; end end % if nargin end function [z, r, residual] = fitcircle(x, varargin) %FITCIRCLE least squares circle fit % % [Z, R] = FITCIRCLE(X) fits a circle to the N points in X minimising % geometric error (sum of squared distances from the points to the fitted % circle) using nonlinear least squares (Gauss Newton) % Input % X : 2xN array of N 2D points, with N >= 3 % Output % Z : center of the fitted circle % R : radius of the fitted circle % % [Z, R] = FITCIRCLE(X, 'linear') fits a circle using linear least % squares minimising the algebraic error (residual from fitting system % of the form ax'x + b'x + c = 0) % % [Z, R] = FITCIRCLE(X, Property, Value, ...) allows parameters to be % passed to the internal Gauss Newton method. Property names can be % supplied as any unambiguous contraction of the property name and are % case insensitive, e.g. FITCIRCLE(X, 't', 1e-4) is equivalent to % FITCIRCLE(X, 'tol', 1e-4). Valid properties are: % % Property: Value: % -------------------------------- % maxits positive integer, default 100 % Sets the maximum number of iterations of the Gauss Newton % method % % tol positive constant, default 1e-5 % Gauss Newton converges when the relative change in the solution % is less than tol % % [X, R, RES] = fitcircle(...) returns the 2 norm of the residual from % the least squares fit % % Example: % x = [1 2 5 7 9 3; 7 6 8 7 5 7]; % % Get linear least squares fit % [zl, rl] = fitcircle(x, 'linear') % % Get true best fit % [z, r] = fitcircle(x) % % Reference: "Least-squares fitting of circles and ellipses", W. Gander, % G. Golub, R. Strebel - BIT Numerical Mathematics, 1994, Springer % This implementation copyright Richard Brown, 2007, but is freely % available to copy, use, or modify as long as this line is maintained error(nargchk(1, 5, nargin, 'struct')) % Default parameters for Gauss Newton minimisation params.maxits = 100; params.tol = 1e-5; % Check x and get user supplied parameters [x, fNonlinear, params] = parseinputs(x, params, varargin{:}); % Convenience variables m = size(x, 2); x1 = x(1, :)'; x2 = x(2, :)'; % 1) Compute best fit w.r.t. algebraic error using linear least squares % % Circle is represented as a matrix quadratic form % ax'x + b'x + c = 0 % Linear least squares estimate found by minimising Bu = 0 s.t. norm(u) = 1 % where u = [a; b; c] % Form the coefficient matrix B = [x1.^2 + x2.^2, x1, x2, ones(m, 1)]; % Least squares estimate is right singular vector corresp. to smallest % singular value of B [U, S, V] = svd(B); u = V(:, 4); % For clarity, set the quadratic form variables a = u(1); b = u(2:3); c = u(4); % Convert to centre/radius z = -b / (2*a); r = sqrt((norm(b)/(2*a))^2 - c/a); % 2) Nonlinear refinement to miminise geometric error, and compute residual if fNonlinear [z, r, residual] = fitcircle_geometric(x, z, r); else residual = norm(B * u); end % END MAIN FUNCTION BODY % NESTED FUNCTIONS function [z, r, residual] = fitcircle_geometric(x, z0, r0) % Use a simple Gauss Newton method to minimize the geometric error fConverged = false; % Set initial u u = [z0; r0]; % Delta is the norm of current step, scaled by the norm of u delta = inf; nIts = 0; for nIts = 1:params.maxits % Find the function and Jacobian [f, J] = sys(u); % Solve for the step and update u h = -J \ f; u = u + h; % Check for convergence delta = norm(h, inf) / norm(u, inf); if delta < params.tol fConverged = true; break end end if ~fConverged warning('fitcircle:FailureToConverge', ... 'Gauss Newton iteration failed to converge'); end z = u(1:2); r = u(3); f = sys(u); residual = norm(f); function [f, J] = sys(u) %SYS Nonlinear system to be minimised - the objective %function is the distance to each point from the fitted circle %contained in u % Objective function f = (sqrt(sum((repmat(u(1:2), 1, m) - x).^2)) - u(3))'; % Jacobian denom = sqrt( (u(1) - x1).^2 + (u(2) - x2).^2 ); J = [(u(1) - x1) ./ denom, (u(2) - x2) ./ denom, repmat(-1, m, 1)]; end % sys end % fitcircle_geometric % END NESTED FUNCTIONS end % fitcircle function [x, fNonlinear, params] = parseinputs(x, params, varargin) % Make sure x is 2xN where N > 3 if size(x, 2) == 2 x = x'; end if size(x, 1) ~= 2 error('fitcircle:InvalidDimension', ... 'Input matrix must be two dimensional') end if size(x, 2) < 3 error('fitcircle:InsufficientPoints', ... 'At least 3 points required to compute fit') end % determine whether we are measuring geometric error (nonlinear), or % algebraic error (linear) fNonlinear = true; switch length(varargin) % No arguments means a nonlinear least squares with defaul parameters case 0 return % One argument can only be 'linear', specifying linear least squares case 1 if strncmpi(varargin{1}, 'linear', length(varargin{1})) fNonlinear = false; return else error('fitcircle:UnknownOption', 'Unknown Option') end % Otherwise we're left with user supplied parameters for Gauss Newton otherwise if rem(length(varargin), 2) ~= 0 error('fitcircle:propertyValueNotPair', ... 'Additional arguments must take the form of Property/Value pairs'); end % Cell array of valid property names properties = {'maxits', 'tol'}; while length(varargin) ~= 0 property = varargin{1}; value = varargin{2}; % If the property has been supplied in a shortened form, lengthen it iProperty = find(strncmpi(property, properties, length(property))); if isempty(iProperty) error('fitcircle:UnkownProperty', 'Unknown Property'); elseif length(iProperty) > 1 error('fitcircle:AmbiguousProperty', ... 'Supplied shortened property name is ambiguous'); end % Expand property to its full name property = properties{iProperty}; switch property case 'maxits' if value <= 0 error('fitcircle:InvalidMaxits', ... 'maxits must be an integer greater than 0') end params.maxits = value; case 'tol' if value <= 0 error('fitcircle:InvalidTol', ... 'tol must be a positive real number') end params.tol = value; end % switch property varargin(1:2) = []; end % while end % switch length(varargin) end %parseinputs function choice = timeoutdlg(delay) prompt='Manual over ride [y/n]?'; f1 = findall(0, 'Type', 'figure'); delay t = timer('TimerFcn', {@closeit f1}, 'StartDelay', delay); start(t); % Call the dialog retvals = inputdlg(prompt,'',1,{'n'}); if numel(retvals) == 0 choice = 'n'; else choice = retvals{1}; end % Delete the timer if strcmp(t.Running, 'on') stop(t); end delete(t); function closeit(src, event, f1) f2 = findall(0, 'Type', 'figure'); fnew = setdiff(f2, f1); if ishandle(fnew); close(fnew); end end end %% function [shape_curvature] = curvature1_v5(Bk,Pcurv1) %%***********************************************************************% %* Curvature measure *% %* Measure shape properties of loops. *% %* *% %* Original author: Dr. Meghan Driscoll, Preetham Manjunatha *% %* Modified by: Bedri Abubaker-Sharif *% %************************************************************************% % % Usage: [shape, Icurv] = curvature(image, boundaryPoint, curvatureThresh, ... % bp_tangent, interp_resolution, loopclose) % Inputs: % Inputs: % binaryImage: input binary image % Bk: for closed curves, this should be a circular list of spaced boundary points % Pcurv1: packed variables for settings; % Outputs: % shape % .curvature - the boundary curvature at each boundary point (uses snakeNum) % Curvatures above or below a cutoff are given the magnitude of the cutoff % .meanNegCurvature - the mean negative curvature % .numIndents - the number of boundary regions over which the curvature is negative % .tangentAngle - the angle of the tangent to the boundary at each boundary point % .tortuosity - the boundary tortuousity (a measure of how bendy the boundary is) % .uncutCurvature - the uncut boundary curvature at each boundary point (uses snakeNum) % Icurv - Output image (padded image to make 3 channel. This fixes the plot) %-------------------------------------------------------------------------- % unpack structure variables v2struct(Pcurv1); % % Perimeter of the binary component % stats = regionprops(binaryimage,'perimeter'); % perimeter = cat(1,stats(1).Perimeter); shape_XY = fliplr(Bk); % convert to (x,y) (col, row) % shape_XY is circular if loopclose == 1 % loop is closed % if there are 5 unique bp: 1 2 3 4 5 % shape_XY: 1 2 3 4 5 1 % bp_positions: 5 1 2 3 4 5 1 % Nb is # of unique boundary points; M is length of shape_XY Nb = size(shape_XY,1) - 1; bp_positions = [shape_XY(Nb,:); shape_XY(1:Nb+1,:)]; % Mth term is already a repeat of 1st boundary point else % loop is open % if there are 5 unique bp: 1 2 3 4 5 % repeat neighboring elements for end points % shape_XY: 1 2 3 4 5 % bp_positions: 2 1 2 3 4 5 4 % Nb is # of unique boundary points; M is length of shape_XY Nb = size(shape_XY,1) ; %bp_positions = [shape_XY(2,:); shape_XY(1:M,:); shape_XY(M-1,:)]; bp_positions = [shape_XY(2,:); shape_XY(1:Nb,:); shape_XY(Nb-1,:)]; end % initialize variables shape_curvature = NaN(Nb,1); shape_uncutCurvature = NaN(Nb,1); shape_meanNegCurvature = NaN; shape_numIndents = NaN; shape_tortuosity = NaN; shape_tangentAngle = NaN(Nb,1); % calculate the curvature (by finding the radius of the osculating circle using three neaby boundary points) for j = 1:Nb % assign the three points that the circle will be fit to such that the slopes are not infinite p1 = bp_positions(j,:); p2 = bp_positions(j+1,:); pstar = p2; p3 = bp_positions(j+2,:); s12 = (p1(2)-p2(2))/(p1(1)-p2(1)); s23 = (p2(2)-p3(2))/(p2(1)-p3(1)); % if exactly one of the two slopes is infinite, then relabel points if isinf(abs(s12)) && ~isinf(abs(s23)) % switch point 2,3 for slope computations point0 = p2; p2 = p3; p3 = point0; s12 = (p1(2)-p2(2))/(p1(1)-p2(1)); % this slope does not change s23 = (p2(2)-p3(2))/(p2(1)-p3(1)); end if isinf(abs(s23)) && ~isinf(abs(s12)) % switch point 1,2 for slope computations point0 = p1; p1 = p2; p2 = point0; % this slope does not change s12 = (p1(2)-p2(2))/(p1(1)-p2(1)); % this slope changes s23 = (p2(2)-p3(2))/(p2(1)-p3(1)); end % if the boundary is flat (s12, s23 both 0, both Inf, both -Inf) if s12==s23 shape_curvature(j) = 0; % if the boundary is curved else % calculate the curvature % mp12: midpoint between line segment 12 % mp23: midpoint between " " 23 % mp13: midpoint between " " 13 mp12 = (p1+p2)/2; mp23 = (p2+p3)/2; mp13 = (p1+p3)/2; % x coordinate of center of circumscribed circle xc = -s23*(mp12(1))/(-s23+s12) + s12*mp23(1)/(-s23+s12) + s12*s23*(mp23(2)-mp12(2))/(-s23+s12); if isequal(s12,0) yc = (-1/s23)*(xc-mp23(1))+mp23(2); else yc = (-1/s12)*(xc-mp12(1))+mp12(2); end shape_uncutCurvature(j) = 1/sqrt((p2(1)-xc)^2+(p2(2)-yc)^2); shape_curvature(j) = shape_uncutCurvature(j); % cutoff the curvature (for visualization) % if shape_curvature(j) > curvatureThresh % shape_curvature(j) = curvatureThresh; % end % determine if the curvature is positive or negative % [In, On] = inpolygon(mp13(1),mp13(2),shape_XY(:,1),shape_XY(:,2)); % % if ~In % shape_curvature(j) = -1*shape_curvature(j); % shape_uncutCurvature(j) = -1*shape_uncutCurvature(j); % end % determine if the curvature is positive or negative % compute tiny line segment from perimeter towards center of fitted circle % find if point lies outside the boundary (curvature is negative) x0 = pstar(1); y0 = pstar(2); xp = x0 + 1e-3*(xc-x0); %order of substraction shows direction yp = y0 + 1e-3*(yc-y0); In = inpolygon(xp,yp,shape_XY(:,1),shape_XY(:,2)); % 1 if xp,yp is inside or on polygon if ~In shape_curvature(j) = -1*shape_curvature(j); shape_uncutCurvature(j) = -1*shape_uncutCurvature(j);; end end end % find the mean negative curvature (really this should use a constant dist snake) listCurve = shape_uncutCurvature(1:Nb); listNegCurve = abs(listCurve(listCurve < 0)); if ~isempty(listNegCurve) shape_meanNegCurvature = sum(listNegCurve)/(Nb); else shape_meanNegCurvature = 0; end % find the number of negative boundary curvature regions curveMask = (listCurve < 0); curveMaskLabeled = bwlabel(curveMask); numIndents = max(curveMaskLabeled); if curveMask(1) && curveMask(end) numIndents = numIndents-1; end shape_numIndents = numIndents; % find the tortuosity (should correct units) shape_tortuosity = sum(gradient(shape_uncutCurvature(1:Nb)).^2)/perimeter; % calculate the direction of the tangent to the boundary if loopclose == 1 %Nb = M-1 %bp_positions_tangent = [shape_XY(M-1,:); shape_XY(1:M,:)]; bp_positions_tangent = [shape_XY(Nb,:); shape_XY(1:Nb+1,:)]; else %Nb = M %bp_positions_tangent = [shape_XY(2,:); shape_XY(1:M,:); shape_XY(M-1,:)]; bp_positions_tangent = [shape_XY(2,:); shape_XY(1:Nb,:); shape_XY(Nb-1,:)]; end for j=1:Nb pt1 = bp_positions_tangent(j,:); pt2 = bp_positions_tangent(j+2,:); shape_tangentAngle(1,j) = mod(atan2(pt1(2)-pt2(2), pt1(1)-pt2(1)), pi); end % shape.XY = shape_XY; % shape.curvature = shape_curvature; % shape.uncutCurvature = shape_uncutCurvature; % shape.meanNegCurvature = shape_meanNegCurvature; % shape.numIndents = shape_numIndents; % shape.tortuosity = shape_tortuosity; % shape.tangentAngle = shape_tangentAngle; end %-------------------------------------------------------------------------- function [xi,yi] = snakeinterp_ba1(x,y,RES) % SNAKEINTERP1 Interpolate the snake to have equal distance RES % [xi,yi] = snakeinterp(x,y,RES) % % RES: resolution desired % update on snakeinterp after finding a bug % Chenyang Xu and Jerry L. Prince, 3/8/96, 6/17/97 % Copyright (c) 1996-97 by Chenyang Xu and Jerry L. Prince % image Analysis and Communications Lab, Johns Hopkins University % convert to column vector % make sure x, y are circular lists (last entry is repeat of first) x = x(:); y = y(:); N = length(x); dx = diff(x); dy = diff(y); d = sqrt(dx.*dx+dy.*dy); % compute the distance from previous node for point 2:N+1 d = [0;d]; % point 1 to point 1 is 0 % now compute the arc length of all the points to point 1 % we use matrix multiply to achieve summing M = length(d); csumd = cumsum(d); % now ready to reparametrize the closed curve in terms of arc length perim = csumd(M); if (perim/RES<3) error('RES too big compare to the length of original curve'); end di = [0:RES:perim]'; xi = interp1(csumd,x,di); yi = interp1(csumd,y,di); N = length(xi); if (perim - di(length(di)) <RES/2) % deal with end boundary condition xi = xi(1:N-1); yi = yi(1:N-1); end % make the list circular xi = [xi;xi(1)]; yi = [yi;yi(1)]; end function C = plotcurvature1_v6(Irgb,boundaries,Bkcurv,Bkpoints,iflag) % raw or SG smoothed boundary points; circular lists Borigs = boundaries; % (row , col ) Xorig = Borigs(:,2); % col indices of original boundary points Yorig = Borigs(:,1); % row indices of original boundary points %% Plot curvature % Bkpoints, Bkcurve are circular lists if iflag == 1 Xspaced = Bkpoints(1:end-1,2); % x coords for spaced boundary points used for curvature computation Yspaced = Bkpoints(1:end-1,1); % y coords for spaced boundary points used for curvature computation Cspaced = Bkcurv(1:end-1); % spaced curvatures % %%% OPTION 1 scattered interpolation to original boundary points F = scatteredInterpolant(Xspaced,Yspaced,Cspaced); C = F(Xorig,Yorig); % circular list of curvatures interpolated to original boundary points else Xspaced = Bkpoints(1:end-1,2); % x coords for spaced boundary points used for curvature computation Yspaced = Bkpoints(1:end-1,1); % y coords for spaced boundary points used for curvature computation Cspaced = Bkcurv(1:end-1); % spaced curvatures %%% OPTION 2 % Bkpoints(double) contains fewer boundary points than Borigs(indices; double if smoothed/interpolated) % Want to interpolate a curve with the same number of boundary points as Borigs % problem: poor X,Y b/c we're only utilizing sparse Bkpoints to reconstruct boundary; % we're not including any info on Borigs, just the number of points in Borigs % so just ignore boundary output (yxinfo) % note: kinfo is circular list; same # of elements as Xorig, Yorig % for circular data: n = # of desired unique points + 1 [~,kinfo] = interpboundary(Bkpoints,size(Borigs,1),Bkcurv); C = kinfo; end % size of data markers in scatter3 dmark = 45; %% overlay interpolated curvature on original and interpolated/sampled boundary points % figure(20); clf; % s1 = subplot(1,2,1); axis equal; hold on % title('k-circum interpolated along raw boundary points') % imshow(Irgb) % scatter3(Xorig,Yorig,zeros(size(Xorig)),[],C,'filled','SizeData',dmark) % % cmap = parula; colormap(cmap); % colorbar; caxis([-0.1 0.1]); % hold off % % s2 = subplot(1,2,2); axis equal; hold on % title('k-circum along smooth,spaced boundary points') % imshow(Irgb) % scatter3(Xspaced,Yspaced,zeros(size(Xspaced)),[],Cspaced,'filled','SizeData',dmark) % % cmap = parula; colormap(cmap); % colorbar; caxis([-0.1 0.1]); % hold off %% % figure(21); clf; % subplot(2,1,1) % plot(C); title('interpolated curvature (k)') % xlabel('k-circum interpolated along raw boundary points') % % subplot(2,1,2) % plot(Cspaced); title('raw curvature (k)') % xlabel('k-circum along smooth,spaced boundary points') % figure; % [x1,y1,z1] = peaks; % surf(x1,y1,z1) end function bdy_sg = sgsmoothboundary(I,bdy,windowWidth,polynomialOrder) %% smooth boundary using savitsky golay filter % make equal spacing before sgsmooth to have each n = size(bdy,1); % circular Borig1 = interpboundary(bdy,n); % circular % get unique points, non-circular x = Borig1(1:end-1,2); y = Borig1(1:end-1,1); % overlap curve of unique boundary points % workaround edge effects problems with sgolayfilt overlap = 2*windowWidth; xplus = [x; x(1:overlap)]; yplus = [y; y(1:overlap)]; % Now smooth with a Savitzky-Golay sliding polynomial filter xs = sgolayfilt(xplus, polynomialOrder, windowWidth); ys = sgolayfilt(yplus, polynomialOrder, windowWidth); %xs = [xs; xs(1)]; ys = [ys; ys(1)]; Borigs1 = [ys,xs]; % keep region with good smoothing, and keep length same as unique boundary points Borigs2 = Borigs1(windowWidth+1:end-windowWidth,:); % circular shift to match indices to original list of boundary points Borigs3 = circshift(Borigs2,windowWidth); % make circular bc currently NONE of the lists are circular bdy_sg = [Borigs3;Borigs3(1,:)]; %Borigs4 = interpboundary(Borigs3c,length(Borig)); Borigs4 = bdy_sg; % f50 = figure(50); %clf(f50); % imshow(I); % axis equal, hold on; % plot(bdy(1,2),bdy(1,1),'r*','LineWidth',10) % plot(bdy(end,2),bdy(end,1),'wx','LineWidth',10) % plot(bdy(:,2),bdy(:,1),'g*','LineWidth',2) % plot(Borigs4(1,2),Borigs4(1,1),'r*','LineWidth',10) % plot(Borigs4(end,2),Borigs4(end,1),'wx','LineWidth',10) % plot(Borigs4(:,2),Borigs4(:,1),'b*','LineWidth',2) % title('Original (green), Savitzky-Golay smoothed (blue)',... % 'FontSize',10); % f50.Position(3:4) = [500 500]; % drawnow % % original boundary % plot(bdy(1,2),bdy(1,1),'r*','LineWidth',10) % plot(bdy(end,2),bdy(end,1),'wx','LineWidth',10) % plot(bdy(:,2),bdy(:,1),'g*','LineWidth',2) % % interpolated boundary for equal spacing % plot(Borig1(1,2),Borig1(1,1),'r*','LineWidth',10) % plot(Borig1(end,2),Borig1(end,1),'wx','LineWidth',10) % plot(Borig1(:,2),Borig1(:,1),'g*','LineWidth',2) % title('Original boundary (green), equal spaced boundary (blue)') % hold off % drawnow end function Bkpoints = defBkpoints(I,Borigs,defBk,Bk_dist,Bk_num) %% boundary point adjustments for curvature estimation % Bk_num = # of UNIQUE boundary points (non-circular) % define spaced boundary points for curvature estimation % Borig is (y,x) circular list or approximately circular (if using sgolay smoothed) % Bkpoints: boundary points (row,col) for curvature computation % assume closed loop: Bkpoints is a circular list if strcmp(defBk,'distance') [xi1,yi1] = snakeinterp_ba1(Borigs(:,2),Borigs(:,1),Bk_dist); Bkpoints = [yi1,xi1]; elseif strcmp(defBk,'number') % will be circular % n = # unique + 1; Bkpoints = interpboundary(Borigs,Bk_num+1); % adjust X,Y coords if last point is too close to first point % only relevant for Bk points where we want equally spaced boundary pfirst = Bkpoints(1,:); psecond = Bkpoints(2,:); plast = Bkpoints(end,:); gaptest = norm(pfirst-plast); gap1 = norm(pfirst-psecond); if gaptest < 1e-8 % Bkpoints is perfectly circular, just finalize Bkpoints(end,:) = Bkpoints(1,:); elseif gaptest < 0.5*gap1 Bkpoints = Bkpoints(1:end-1,:); Bkpoints = [Bkpoints; Bkpoints(1,:)]; end else Bkpoints = Borigs; end % f10 = figure(10); clf(f10); % imshow(I); % axis equal, hold on % plot(Bkpoints(1,2),Bkpoints(1,1),'r*','LineWidth',10) % plot(Bkpoints(:,2),Bkpoints(:,1),'b*','LineWidth',3) % title('Smoothed and spaced boundary') % hold off % drawnow end function varargout = mavgsmoothboundary(mfilter,varargin) %% % Moving average; Should be an odd number s = mfilter; n = length(varargin); % # boundary quantities varargout = cell(1,n); if s>0 % then do movavg for k = 1:n % each boundary quantity bdy = varargin{k}; sbdy = bdy; for i=1:s % sbdy = sbdy+circshift(bdy,i)+circshift(bdy,-i); end sbdy = sbdy/(2*s+1); varargout{k} = sbdy; end % figure() %plot(bdy(1:200,1),bdy(1:200,2),'x',sbdy(1:200,1),sbdy(1:200,2),'.r') else varargout = varargin; end end function K = curve_fitcircle(bdy,delta) % INPUT % bdy: list of (y,x) or (row,col) positions of boundary points (circular list) % delta is # of points to the right/left to consider fitting % OUTPUT % C1: local curvatures based on fitting a circle that minimizes geometric error % same length as bdy (circular list) % make equal spaced boundary segments for consistency to delta around the perimeter n = size(bdy,1); % keep same length as original boundary bdyI = interpboundary(bdy,n); % circular B = fliplr(bdyI); % now it's [x, y] coords Pn = B(1:end-1,:); % only uniqe boundary points considered for circshift to work properly curve = zeros(length(Pn),1); for i=1:length(Pn) s=circshift(Pn,delta-i+1);%not delta-i % Fit to a circle the 2*delta+1 points centered at perimeter index "i" [zz,rr]=fitcircle(s(1:2*delta+1,:)'); % not -1 xc = zz(1); yc = zz(2); % circle center coords (x, y) x0 = Pn(i,1); y0 = Pn(i,2); % boundary index coords (x, y) if isnan(rr) % rr is NaN; % so curvature is 1/Inf = zero curve(i)=0; else curve(i) = 1/rr; % unsigned curvature % determine if the curvature is positive or negative % compute tiny line segment from perimeter towards center of fitted circle % find if point lies outside the boundary (curvature is negative) xp = x0 + 1e-3*(xc-x0); %order of substraction shows direction yp = y0 + 1e-3*(yc-y0); In = inpolygon(xp,yp,Pn(:,1),Pn(:,2)); % 1 if xp,yp is inside or on polygon if ~In curve(i) = -curve(i); end end end % make circular K = [curve; curve(1,:)]; %ncurve=curve/(abs(mean(curve))); end %--------------- function C2 = plotcurvature2_v1(Irgb,boundaries,c2a,Bkpoints,iflag) % varargin when using spaced points: Bkpoints, iflag % otherwise only 3 inputs % raw or SG smoothed boundary points; circular lists Borigs = boundaries; % (row , col ) Xorig = Borigs(:,2); % col indices of original boundary points Yorig = Borigs(:,1); % row indices of original boundary points %% Plot curvature % Bkpoints, Bkcurve are circular lists if nargin > 3 if iflag == 1 Xspaced = Bkpoints(1:end-1,2); % x coords for spaced boundary points used for curvature computation Yspaced = Bkpoints(1:end-1,1); % y coords for spaced boundary points used for curvature computation C2spaced = c2a(1:end-1); % spaced curvatures % %%% OPTION 1 scattered interpolation to original boundary points F = scatteredInterpolant(Xspaced,Yspaced,C2spaced); C2 = F(Xorig,Yorig); % circular list of curvatures interpolated to original boundary points else Xspaced = Bkpoints(1:end-1,2); % x coords for spaced boundary points used for curvature computation Yspaced = Bkpoints(1:end-1,1); % y coords for spaced boundary points used for curvature computation C2spaced = c2a(1:end-1); % spaced curvatures %%% OPTION 2 % Bkpoints(double) contains fewer boundary points than Borigs(indices; double if smoothed/interpolated) % Want to interpolate a curve with the same number of boundary points as Borigs % problem: poor X,Y b/c we're only utilizing sparse Bkpoints to reconstruct boundary; % we're not including any info on Borigs, just the number of points in Borigs % so just ignore boundary output (yxinfo) % note: kinfo is circular list; same # of elements as Xorig, Yorig % for circular data: n = # of desired unique points + 1 [~,kinfo] = interpboundary(Bkpoints,size(Borigs,1),c2a); C2 = kinfo; end else C2 = c2a; % stay circular end % size of data markers in scatter3 dmark = 45; %% overlay interpolated curvature on original and interpolated/sampled boundary points % figure(40); % s1 = tiledlayout('flow'); s1.Padding='compact';s1.TileSpacing='compact'; % nexttile(1,[1 3]); axis equal; hold on % if nargin>3 % title('k-fitcircle interpolated along spaced boundary points') % else % title('k-fitcircle along boundary points') % end % imshow(Irgb) % f40ax1 = scatter3(Xorig,Yorig,zeros(size(Xorig)),[],C2,'filled','SizeData',dmark); % cmap = parula; colormap(cmap); % colorbar; f40ax1.Parent.CLim = [-0.22 0.22]; %caxis([-0.03 0.03]); % hold off % % if nargin>3 % nexttile(4,[1 3]); axis equal; hold on % title('k-fitcircle along spaced boundary points') % imshow(Irgb) % f40ax2 = scatter3(Xspaced,Yspaced,zeros(size(Xspaced)),[],C2spaced,'filled','SizeData',dmark); % cmap = parula; colormap(cmap); % colorbar; f40ax2.Parent.CLim = [-0.03, 0.03]; %caxis([-0.03 0.03]); % hold off % end %% % figure(41); % s2 = tiledlayout('flow'); s2.Padding='compact';s2.TileSpacing='compact'; % nexttile(1,[1 3]); % plot(C2); % if nargin > 3 % title('interpolated curvature (k)') % xlabel('k-fitcircle interpolated along raw boundary points') % else % title('curvature (k-fitcircle)') % xlabel('k-fitcircle along raw boundary points') % end % % if nargin>3 % nexttile(4, [1 3]); % plot(C2spaced); title('raw curvature (k)') % xlabel('k-fitcircle along smooth,spaced boundary points') % end % figure; % [x1,y1,z1] = peaks; % surf(x1,y1,z1) end function [A , c] = MinVolEllipse(P, tolerance) % [A , c] = MinVolEllipse(P, tolerance) % Finds the minimum volume enclsing ellipsoid (MVEE) of a set of data % points stored in matrix P. The following optimization problem is solved: % % minimize log(det(A)) % subject to (P_i - c)' * A * (P_i - c) <= 1 % % in variables A and c, where P_i is the i-th column of the matrix P. % The solver is based on Khachiyan Algorithm, and the final solution % is different from the optimal value by the pre-spesified amount of 'tolerance'. % % inputs: %--------- % P : (d x N) dimnesional matrix containing N points in R^d. % tolerance : error in the solution with respect to the optimal value. % % outputs: %--------- % A : (d x d) matrix of the ellipse equation in the 'center form': % (x-c)' * A * (x-c) = 1 % c : 'd' dimensional vector as the center of the ellipse. % % example: % -------- % P = rand(5,100); % [A, c] = MinVolEllipse(P, .01) % % To reduce the computation time, work with the boundary points only: % % K = convhulln(P'); % K = unique(K(:)); % Q = P(:,K); % [A, c] = MinVolEllipse(Q, .01) % % % Nima Moshtagh (nima@seas.upenn.edu) % University of Pennsylvania % % December 2005 % UPDATE: Jan 2009 %%%%%%%%%%%%%%%%%%%%% Solving the Dual problem%%%%%%%%%%%%%%%%%%%%%%%%%%%5 % --------------------------------- % data points % ----------------------------------- [d N] = size(P); Q = zeros(d+1,N); Q(1:d,:) = P(1:d,1:N); Q(d+1,:) = ones(1,N); % initializations % ----------------------------------- count = 1; err = 1; u = (1/N) * ones(N,1); % 1st iteration % Khachiyan Algorithm % ----------------------------------- while err > tolerance, X = Q * diag(u) * Q'; % X = \sum_i ( u_i * q_i * q_i') is a (d+1)x(d+1) matrix M = diag(Q' * inv(X) * Q); % M the diagonal vector of an NxN matrix [maximum j] = max(M); step_size = (maximum - d -1)/((d+1)*(maximum-1)); new_u = (1 - step_size)*u ; new_u(j) = new_u(j) + step_size; count = count + 1; err = norm(new_u - u); u = new_u; end %%%%%%%%%%%%%%%%%%% Computing the Ellipse parameters%%%%%%%%%%%%%%%%%%%%%% % Finds the ellipse equation in the 'center form': % (x-c)' * A * (x-c) = 1 % It computes a dxd matrix 'A' and a d dimensional vector 'c' as the center % of the ellipse. U = diag(u); % the A matrix for the ellipse % -------------------------------------------- A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' ); % center of the ellipse % -------------------------------------------- c = P * u; end function Ellipse_plot(A, C, varargin) % % Ellipse_Plot(A,C,N) plots a 2D ellipse or a 3D ellipsoid % represented in the "center" form: % % (x-C)' A (x-C) <= 1 % % A and C could be the outputs of the function: "MinVolEllipse.m", % which computes the minimum volume enclosing ellipsoid containing a % set of points in space. % % Inputs: % A: a 2x2 or 3x3 matrix. % C: a 2D or a 3D vector which represents the center of the ellipsoid. % N: the number of grid points for plotting the ellipse; Default: N = 20. % % Example: % % P = rand(3,100); % t = 0.001; % [A , C] = MinVolEllipse(P, t) % figure % plot3(P(1,:),P(2,:),P(3,:),'*') % hold on % Ellipse_plot(A,C) % % % Nima Moshtagh % nima@seas.upenn.edu % University of Pennsylvania % Feb 1, 2007 % Updated: Feb 3, 2007 %%%%%%%%%%% start %%%%%%%%%%%%%%%%%%%%%%%%%%% N = 20; % Default value for grid % See if the user wants a different value for N. %---------------------------------------------- if nargin > 2 N = varargin{1}; end % check the dimension of the inputs: 2D or 3D %-------------------------------------------- if length(C) == 3, Type = '3D'; elseif length(C) == 2, Type = '2D'; else display('Cannot plot an ellipse with more than 3 dimensions!!'); return end % "singular value decomposition" to extract the orientation and the % axes of the ellipsoid [U D V] = svd(A); if strcmp(Type, '2D'), % get the major and minor axes %------------------------------------ a = 1/sqrt(D(1,1)); b = 1/sqrt(D(2,2)); theta = [0:1/N:2*pi+1/N]; % Parametric equation of the ellipse %---------------------------------------- state(1,:) = a*cos(theta); state(2,:) = b*sin(theta); % Coordinate transform %---------------------------------------- X = V * state; X(1,:) = X(1,:) + C(1); X(2,:) = X(2,:) + C(2); elseif strcmp(Type,'3D'), % generate the ellipsoid at (0,0,0) %---------------------------------- a = 1/sqrt(D(1,1)); b = 1/sqrt(D(2,2)); c = 1/sqrt(D(3,3)); [X,Y,Z] = ellipsoid(0,0,0,a,b,c,N); % rotate and center the ellipsoid to the actual center point %------------------------------------------------------------ XX = zeros(N+1,N+1); YY = zeros(N+1,N+1); ZZ = zeros(N+1,N+1); for k = 1:length(X), for j = 1:length(X), point = [X(k,j) Y(k,j) Z(k,j)]'; P = V * point; XX(k,j) = P(1)+C(1); YY(k,j) = P(2)+C(2); ZZ(k,j) = P(3)+C(3); end end end % Plot the ellipse %---------------------------------------- if strcmp(Type,'2D'), plot(X(1,:),X(2,:),'c','LineWidth',2); hold on; plot(C(1),C(2),'m*','MarkerSize',10,LineWidth=3); axis equal, grid else mesh(XX,YY,ZZ); axis equal hidden off end end function [aspratio,A0,A0c,a0majax] = aspectratio(bdyraw) % prm0 = v2struct(minfont,movieFrame,do_plotEllipseries) %v2struct(prm0); %byraw = rows - bdyraw(:,1) + 1; byraw = bdyraw(:,1); bxraw = bdyraw(:,2); xpts = [bxraw,byraw]; %xmean = xpts-mean(xpts); % if do_plotEllipseries % f2 = figure(2); imshow(Irgb); axis equal; hold on; % plot(xpts(:,1),xpts(:,2),'r','LineWidth',3); % FontSize = minfont; %min(minfont,floor(min(cols,rows)/10)); % text(6,6,num2str(movieFrame),'Color','white','FontSize',FontSize) % % set(gca,'YDir','reverse'); % % [A0,A0c] = MinVolEllipse(xpts',0.01); % a0 = sqrt(1/min(A0(1),A0(4))); % b0 = sqrt(1/max(A0(1),A0(4))); % eccent = sqrt(1-(b0/a0)^2); % aspratio = a0/b0; % %stats.Eccentricity % % Ellipse_plot(A0,A0c); % varargout{1} = f2; % else [A0,A0c] = MinVolEllipse(xpts',0.01); a0 = sqrt(1/min(A0(1),A0(4))); b0 = sqrt(1/max(A0(1),A0(4))); eccent = sqrt(1-(b0/a0)^2); aspratio = a0/b0; a0majax = 2*a0; % % end end function v = recordvid(currentframe,vidparams) % save captured frame as avi % recall vidparams = v2struct(fname,vidrate,basystem) v2struct(vidparams) if strcmp('linux',basystem) v = VideoWriter(fname,'Motion JPEG AVI'); else v = VideoWriter(fname,'MPEG-4'); end v.FrameRate = vidrate; %playback framerate open(v); axis equal writeVideo(v,currentframe.cdata); %pause(0.1) close(v); end