diff --git a/BasicPlotMaps.m b/BasicPlotMaps.m index 9196801..f8e5c37 100644 --- a/BasicPlotMaps.m +++ b/BasicPlotMaps.m @@ -18,6 +18,15 @@ function BasicDrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,is barratio = 0.02; numK = double(max(double(numK),double(max(gIX)))); +% down-sample +displaymax = 1000; +numcell = size(M,1); +if numcell > displaymax, + skip = round(numcell/displaymax); + M = M(1:skip:end,:); + gIX = gIX(1:skip:end,:); +end + % sort traces by index im = M; [nLines,nFrames] = size(im); diff --git a/DrawClusters.m b/DrawClusters.m index 964da88..524ce77 100644 --- a/DrawClusters.m +++ b/DrawClusters.m @@ -1,201 +1,293 @@ -function DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite) -pos = get(gca,'Position'); -barratio = 0.03; +function DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite,isPopout,isPlotLines,isPlotFictive) +axis off; +pos = get(h1,'Position'); % [left, bottom, width, height] -%% Prepare cluster data -% down-sample -displaymax = 1000; -numcell = size(M,1); -if numcell > displaymax, - skip = round(numcell/displaymax); - M = M(1:skip:end,:); - gIX = gIX(1:skip:end,:); -end - -numK = double(max(double(numK),double(max(gIX)))); +if isPlotLines, + if 1, % (just to collapse isPlotLines = true) + %% settings + len = pos(3); + fps = 1.97; + + % set position grid + nLines = length(unique(gIX)); + if isPlotFictive, + nExtraRows = 3; % number of extra rows + else + nExtraRows = 2; + end + nRows = max(8,nLines)+nExtraRows; + lineH = pos(4)/nRows; +% Lpos = [pos(1),pos(2)+pos(4)-lineH*nRows,pos(3),lineH*nRows]; +% % Lpos = [pos(1),pos(2)+pos(4)-lineH*(nLines+nExtraRows),pos(3),lineH*(nLines+nExtraRows)]; +% set(h1,'Position',Lpos); -% sort traces by index -im = M; -[nLines,nFrames] = size(im); -[~,I] = sort(gIX); -im = im(I,:); + % figure('Position',[500,900-(min(nLines,15)+nDiv)*50,600,(min(nLines,15)+nDiv)*50],'color',[1 1 1]); + xv = 1:size(M,2); + + %% Set colormap + if strcmp(clrmap,'jet'), + clr = flipud(jet(double(numK)))*0.8; % make darker by *0.8 + else % if strcmp(clrmap,'hsv'), + clr = hsv(round(double(numK)*1.1))*0.7; % make darker by *0.7 + end + % clr = lines(nLines); + % clr = ones(nLines,3)*0.3; % uniform charcoal + + %% draw stimulus bar + stimbar = GetStimBar(1,stim); + rpos = [pos(1),pos(2)+pos(4)-3/4*lineH,len,lineH/2]; + h = axes('Position',rpos); + + image(stimbar); + set(h,'Xtick',[],'Ytick',[]) + + % draw fish icon + rpos = [pos(1)*0.2,pos(2)+pos(4)-3/4*lineH,pos(1)*0.7,lineH/2]; + axes('Position',rpos); + DrawFishIcon; + + %% draw cluster means + Ymeans = zeros(nLines,size(M,2)); + U = unique(gIX); + for j = 1:nLines, + k = U(j); + % ymean = Cz(k,:); + Ys = M(gIX==k,:); + ymean = mean(Ys,1); + Ymeans(k,:) = ymean; + rpos = [pos(1),pos(2)+pos(4)-lineH*(j+1),len,lineH]; + axes('Position',rpos); hold on; + + ySTD=std(Ys); % /sqrt(size(Ys,1)) + ySTD_upper=ymean+ySTD; + ySTD_lower=ymean-ySTD; + h = fill([xv fliplr(xv)], [ySTD_upper fliplr(ySTD_lower)],0.5+0.3*clr(k,:)); + set(h, 'EdgeColor','none') + + plot(xv,ymean,'-','Linewidth',1,'color',clr(k,:)) + axis tight;axis off + end + + %% plot scale bar + axes('Position',[pos(1),pos(2)+pos(4)-lineH*(nLines+2),len,lineH]); hold on; + if size(M,2)<1200, % <~10 min, plot 1min scale bar + plot([1,60*fps],[0.75,0.75],'k-'); + text(pos(1)+60*fps/2,0.5,'1 min','HorizontalAlignment','center') + else % >~10 min, plot 10min scale bar + plot([1,600*fps],[0.75,0.75],'k-'); + text(pos(1)+600*fps/2,0.5,'10 min','HorizontalAlignment','center') + end + xlim([xv(1),xv(end)]);ylim([0,1]); + axis off + + %% draw fictive behavior bar + if isPlotFictive, + h = axes('Position',[pos(1),pos(2)+pos(4)-lineH*(nLines+3),len,0.7/(nLines+nExtraRows)]); + DrawFictiveBar(h,fictive); + + axes('Position',[0.02,pos(2)+pos(4)-lineH*(nLines+3),pos(1)-0.02,0.7/(nLines+nExtraRows)]); + DrawArrowsIcon(isPopout); + end -%% convert imagesc effect into rgb matrix 'RGB' -if dataFR, % is drawing cell-response fluorescence - cmap = gray(64); - im = AutoScaleImage0to1(im); % scaling to min/max 0/1 - minlim = 0; - maxlim = 1; -else % create blue-white-red map, for regression results - cmap = zeros(64,3); - cmap(:,1) = [linspace(0,1,32), linspace(1,1,32)]; - cmap(:,2) = [linspace(0,1,32), linspace(1,0,32)]; - cmap(:,3) = [linspace(1,1,32), linspace(1,0,32)]; - minlim = -1; %min(min(im)); - maxlim = 1; %max(max(im)); -end -RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap + end + +else -%% add vertical color code -if exist('clrmap','var'), - if strcmp(clrmap,'jet'), - temp = flipud(jet(numK)); + if ~isPopout, + barratio = 0.025; + else + barratio = 0.05; + end + %% Prepare cluster data + % down-sample + displaymax = 1000; + numcell = size(M,1); + if numcell > displaymax, + skip = round(numcell/displaymax); + M = M(1:skip:end,:); + gIX = gIX(1:skip:end,:); + end + +% numK = double(max(double(numK),double(max(gIX)))); + + % sort traces by index + im = M; + [nLines,nFrames] = size(im); + [~,I] = sort(gIX); + im = im(I,:); + + %% convert imagesc effect into rgb matrix 'RGB' + if dataFR, % is drawing cell-response fluorescence + cmap = gray(64); + im = AutoScaleImage0to1(im); % scaling to min/max 0/1 + minlim = 0; + maxlim = 1; + else % create blue-white-red map, for regression results + cmap = zeros(64,3); + cmap(:,1) = [linspace(0,1,32), linspace(1,1,32)]; + cmap(:,2) = [linspace(0,1,32), linspace(1,0,32)]; + cmap(:,3) = [linspace(1,1,32), linspace(1,0,32)]; + minlim = -1; %min(min(im)); + maxlim = 1; %max(max(im)); + end + RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap + + %% add vertical color code + if exist('clrmap','var'), + if strcmp(clrmap,'jet'), + temp = flipud(jet(numK)); + else % 'hsv' + temp = hsv(round(numK*1.1)); + end else % 'hsv' temp = hsv(round(numK*1.1)); end -else % 'hsv' - temp = hsv(round(numK*1.1)); -end -cmap2 = vertcat(temp(1:numK,:),[0,0,0]); % extend colormap to include black -bwidth = max(round(nFrames/30),1); - -idx = gIX(I); -ix_div = [find(diff(idx));length(idx)]; - -bars = ones(nLines,bwidth); -if numK>1, - bars(1:ix_div(1),:) = idx(ix_div(1)); - for i = 2:length(ix_div), - % paint color - bars(ix_div(i-1)+1:ix_div(i),:) = idx(ix_div(i)); - end -end -im_bars = reshape(cmap2(bars,:),[size(bars),3]); -% add a white division margin -div = ones(nLines,round(bwidth/2),3); -% put 3 parts together -im = horzcat(RGB,div,im_bars); - -%% plot figure - -[s1,s2,s3] = size(im); - -hold off; -set(h1,'Position',[pos(1),pos(2)+pos(4)*barratio*3,pos(3),pos(4)*(1-4*barratio)]); -if s1<30, - temp = im; - im = ones(30,s2,s3); - im(1:s1,:,:) = temp; -end -image(im); -set(gca, 'box', 'off') - -hold on; - -% plot cluster division lines -plot([0.5,s2+0.5],[0.5,0.5],'k','Linewidth',0.5); -if numK>1, - for i = 1:length(ix_div),% = numK-1, - y = ix_div(i)+0.5; - plot([0.5,s2+0.5],[y,y],'k','Linewidth',0.5); + cmap2 = vertcat(temp(1:numK,:),[0,0,0]); % extend colormap to include black + bwidth = max(round(nFrames/30),1); + + idx = gIX(I); + ix_div = [find(diff(idx));length(idx)]; + + bars = ones(nLines,bwidth); + if numK>1, + bars(1:ix_div(1),:) = idx(ix_div(1)); + for i = 2:length(ix_div), + % paint color + bars(ix_div(i-1)+1:ix_div(i),:) = idx(ix_div(i)); + end end -end - -ylabel(['Number of cells: ' num2str(numcell)]); -set(gca,'YTick',[],'XTick',[]); -set(gcf,'color',[1 1 1]); - -colormap gray; - -% write in number label -x = s2*1.003; -y_last = 0; -for i = 1:length(ix_div), - % avoid label crowding - margin = 0.015*s1; - y0 = ix_div(i)+0.5; - y = max(y_last+margin,y0); - if i1, + for i = 1:length(ix_div),% = numK-1, + y = ix_div(i)+0.5; + plot([0.5,s2+0.5],[y,y],'k','Linewidth',0.5); + end + end + + ylabel(['Number of cells: ' num2str(numcell)]); + set(gca,'YTick',[],'XTick',[]); + set(gcf,'color',[1 1 1]); + + colormap gray; + + % write in number label + x = s2*1.003; + y_last = 0; + for i = 1:length(ix_div), + % avoid label crowding + margin = 0.015*s1; + y0 = ix_div(i)+0.5; + y = max(y_last+margin,y0); + if i0 & (cinds+circle_inds)<=dimv_yxz(1)*dimv_yxz(2)); + zinds=dimv_yxz(1)*dimv_yxz(2)*(CInfo(cIX(j)).slice-1); + ix = find(U==gIX(j)); + ixs = cinds+circle_inds(labelinds)+zinds; + ave_stack2(ixs)=cmap(ix,1)*weight + ave_stack2(ixs)*(1-weight); + ixs = cinds+circle_inds(labelinds)+zinds+stacklen; + ave_stack2(ixs)=cmap(ix,2)*weight + ave_stack2(ixs)*(1-weight); + ixs = cinds+circle_inds(labelinds)+zinds+stacklen*2; + ave_stack2(ixs)=cmap(ix,3)*weight + ave_stack2(ixs)*(1-weight); +end + +%% view stack sequentially +% figure; +% for i_plane = 1:nPlanes, +% im = squeeze(ave_stack2(:,:,i_plane,:)); +% image(im); +% axis image; axis off +% pause(0.1) +% end + +% lay out + +h = figure('Position',[10,20,round(1136*1.4),round(683*1.4)],'color','k'); +hold on; + +axes('Position',[0,0,1,1]); % BS fix to get the figure background to stay when saving pics +image(zeros(3,3,3));axis off + +numRow = 3; +numP = nPlanes - mod(nPlanes,numRow); +numCol = numP/numRow; + +for i_plane = 1:numP, + im = squeeze(ave_stack2(:,:,i_plane,:)); + im = imresize(im,0.25); + + [col, row] = ind2sub([numCol,numRow],i_plane); % this is intensionally inverted... + posVec = [(col-1)*1/numCol,1-row*1/numRow,1/numCol,1/numRow]; + axes('Position',posVec); + image(im); + axis image; axis off + % pause(0.1) +end + +end \ No newline at end of file diff --git a/GUI_FishExplorer.m b/GUI_FishExplorer.m index a417875..04405e0 100644 --- a/GUI_FishExplorer.m +++ b/GUI_FishExplorer.m @@ -21,9 +21,10 @@ % global VAR; % load('VAR_current.mat','VAR'); % load('CONSTs_current.mat','CONSTs'); -% hfig = GUI_FishExplorer(CONSTs,VAR); +% data_dir = [pwd '\example data']; +% [hfig,fcns] = GUI_FishExplorer(data_dir,CONSTs,VAR); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Full datasets are stored as 'CONST_F?.mat' (? = fish number), to load from the GUI +% Full datasets are stored as 'CONST_F?_fast.mat' (? = fish number), to load from the GUI %% User Interface: % function hfig = GUI_FishExplorer(CONSTs,VAR,CONST,i_fish) @@ -57,7 +58,11 @@ %% Pass external variables into appdata (stored with main figure handle) setappdata(hfig,'CONSTs',CONSTs); setappdata(hfig,'VAR',VAR); -nFish = length(CONSTs); +nFish = length(CONSTs) + 1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% fish protocol sets (different sets have different parameters) +M_fish_set = [1, 1, 1, 1, 1, 1, 1, 2, 3, 3]; % M = Matrix +setappdata(hfig,'M_fish_set',M_fish_set); % parameters / constants setappdata(hfig,'z_res',19.7); % resoltion in z, um per slice @@ -81,9 +86,6 @@ QuickUpdateFish(hfig,i_fish,'init'); %% Initialize internal params into appdata -% fish protocol sets (different sets have different parameters) -M_fish_set = [1, 1, 1, 1, 1, 1, 1, 2, 1]; % M = Matrix -setappdata(hfig,'M_fish_set',M_fish_set); % thresholds thres_merge = 0.9; @@ -100,6 +102,8 @@ setappdata(hfig,'clrmap','hsv'); setappdata(hfig,'opID',0); setappdata(hfig,'rankID',0); +setappdata(hfig,'isPlotLines',0); +setappdata(hfig,'isPlotFictive',1); setappdata(hfig,'rankscore',[]); setappdata(hfig,'isCentroid',0); setappdata(hfig,'isWkmeans',1); % in autoclustering, with/without kmeans @@ -148,7 +152,7 @@ % various UI element handles ('global' easier than passing around..) global hback hfwd hclassname hclassmenu hclusgroupmenu hclusmenu hclusname... - hdatamenu hopID hcentroid hdataFR hloadfish hfishnum hstimreg hmotorreg... + hdatamenu hopID hdataFR hloadfish hfishnum hstimreg hmotorreg... hcentroidreg; %% UI ----- tab one ----- (General) @@ -176,12 +180,30 @@ 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... 'Callback',@pushbutton_writeZstack_Callback); +i=i+n; +n=2; % plots selected cells on anatomy z-stack, tiled display +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Tile Zstack',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_tileZstack_Callback); + i=i+n; n=2; % make new figure without the GUI components, can save manually from default menu uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Popup plot',... 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... 'Callback',@pushbutton_popupplot_Callback); +i=i+n; +n=2; % popupplot option: whether to plot cluster mean lines instead of all raw traces +uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Plot lines',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@checkbox_isPlotLines_Callback); + +i=i+n; +n=2; % popupplot option: whether to plot fictive behavior bar +uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Plot fictive','Value',1,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@checkbox_isPlotFictive_Callback); + i=i+n; n=3; % export main working variables to workspace, can customize! uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Export to workspace',... @@ -233,7 +255,7 @@ i=i+n; n=4; % only centroids (~mean) of clusters shown on left-side plot, the rest is unchanged -hcentroid = uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Display centroids of clusters',... +uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Display centroids of clusters',... 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... 'Callback',@checkbox_showcentroids_Callback); @@ -413,7 +435,7 @@ i=i+n; n=2; % (unlike stim regressors, names hardcoded, not importet from regressor...) -menu = {'(choose)','right swims','left swims','forward swims','raw right','raw left','raw average'}; +menu = {'(choose)','left swims','right swims','forward swims','raw left','raw right','raw average'}; hmotorreg = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',menu,'Value',1,... 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... 'Callback',@popup_getmotorreg_Callback); @@ -493,19 +515,19 @@ 'Callback',{@pushbutton_IterCentroidRegression_Callback}); %% UI row 3: (TBD) -i_row = 3; % old idea: display correlation values of multiple regressors (instead of fluo. trace) -i = 1;n = 0; % and then can cluster that and further manipulate... - -i=i+n; -n=4; % did not implement combination of regressors in this version... but display option is coded (dataFR==0) -uicontrol('Parent',tab{i_tab},'Style','text','String','regressor combos??(TBD)',... - 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); - -i=i+n; -n=4; % also data storage is coded, 'M' could be loaded from M_0_fluo or M_0_reg (storing reg corr values) -hdataFR = uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Display fluo./regression',... - 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... - 'Callback',@checkbox_dataFR_Callback); +% i_row = 3; % old idea: display correlation values of multiple regressors (instead of fluo. trace) +% i = 1;n = 0; % and then can cluster that and further manipulate... +% +% i=i+n; +% n=4; % did not implement combination of regressors in this version... but display option is coded (dataFR==0) +% uicontrol('Parent',tab{i_tab},'Style','text','String','regressor combos??(TBD)',... +% 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); +% +% i=i+n; +% n=4; % also data storage is coded, 'M' could be loaded from M_0_fluo or M_0_reg (storing reg corr values) +% hdataFR = uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Display fluo./regression',... +% 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... +% 'Callback',@checkbox_dataFR_Callback); %% UI ----- tab four ----- (Clustering etc.) i_tab = 4; @@ -977,6 +999,17 @@ function pushbutton_writeZstack_Callback(hObject,~) close(h) end +function pushbutton_tileZstack_Callback(hObject,~) +hfig = getParentFigure(hObject); +isfullfish = getappdata(hfig,'isfullfish'); +if isfullfish, + disp('Rendering...') + DrawTiledPics(hfig); +else + errordlg('Load full fish first!'); +end +end + function out = makeDisk2(radius, dim) center=floor(dim/2)+1; out=zeros(dim); @@ -997,7 +1030,7 @@ function pushbutton_popupplot_Callback(hObject,~) % figure('Position',[50,100,853,512]); % right plot size: 2048x1706 % h1 = axes('Position',[0, 0, 0.5, 1]); % left ~subplot % h2 = axes('Position',[0.50, 0, 0.5, 1]); % right ~subplot -figure('Position',[50,100,1600,900]); +figure('Position',[50,100,1700,900],'color',[1 1 1]); h1 = axes('Position',[0.05, 0.03, 0.53, 0.94]); % left ~subplot h2 = axes('Position',[0.61, 0.03, 0.35, 0.94]); % right ~subplot @@ -1017,22 +1050,39 @@ function pushbutton_popupplot_Callback(hObject,~) clrmap = getappdata(hfig,'clrmap'); rankscore = getappdata(hfig,'rankscore'); rankID = getappdata(hfig,'rankID'); - iswrite = (rankID>=2); +isPlotLines = getappdata(hfig,'isPlotLines'); +isPlotFictive = getappdata(hfig,'isPlotFictive'); + +isPopout = 1; % left subplot axes(h1); if isCentroid, [C,~] = FindCentroid(gIX,M); - DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); + DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,... + iswrite,isPopout,isPlotLines,isPlotFictive); +% DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,iswrite,ispopout); else - DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); + DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,... + iswrite,isPopout,isPlotLines,isPlotFictive); +% DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite,ispopout); end % right subplot axes(h2); -DrawClustersOnMap_LSh(CInfo,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,clrmap); +DrawClustersOnMap_LSh(CInfo,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,clrmap,'full'); + +end + +function checkbox_isPlotLines_Callback(hObject,~) +hfig = getParentFigure(hObject); +setappdata(hfig,'isPlotLines',get(hObject,'Value')); +end +function checkbox_isPlotFictive_Callback(hObject,~) +hfig = getParentFigure(hObject); +setappdata(hfig,'isPlotFictive',get(hObject,'Value')); end function pushbutton_exporttoworkspace_Callback(hObject,~) @@ -1069,12 +1119,12 @@ function QuickUpdateFish(hfig,new_i_fish,init) %#ok numcell = CONSTs{new_i_fish}.numcell; temp = CONSTs{new_i_fish}.CRAZ; -CRAZ = zeros(numcell,size(temp,2)); -CRAZ(CIX,:) = temp; +CellRespAvr = zeros(numcell,size(temp,2)); +CellRespAvr(CIX,:) = temp; temp = CONSTs{new_i_fish}.CRZt; -CRZt = zeros(numcell,size(temp,2)); -CRZt(CIX,:) = temp; +CellResp = zeros(numcell,size(temp,2)); +CellResp(CIX,:) = temp; temp = CONSTs{new_i_fish}.CInfo; CInfo(numcell).center = ''; @@ -1085,17 +1135,21 @@ function QuickUpdateFish(hfig,new_i_fish,init) %#ok CInfo(numcell).x_minmax = ''; CInfo(CIX) = temp; -setappdata(hfig,'CRAZ',CRAZ); % Matrix array, misc collection -setappdata(hfig,'CRZt',CRZt); % Matrix array, misc collection -setappdata(hfig,'CInfo',CInfo); % Matrix array, misc collection +setappdata(hfig,'CellRespAvr',CellRespAvr); +setappdata(hfig,'CellResp',CellResp); +setappdata(hfig,'CInfo',CInfo); UpdateFishData(hfig,new_i_fish); if ~exist('init','var'), - UpdateFishDisplay(hfig,new_i_fish); + UpdateFishDisplay(hfig); end end -function UpdateFishData(hfig,new_i_fish) +function UpdateFishData(hfig,new_i_fish) % loading steps shared by both quick-load and full-load +M_fish_set = getappdata(hfig,'M_fish_set'); +fishset = M_fish_set(new_i_fish); +setappdata(hfig,'fishset',fishset); + UpdateDatamenu_Direct(hfig,1); % save ClusGroup before updating, if applicable @@ -1122,12 +1176,12 @@ function UpdateFishData(hfig,new_i_fish) setappdata(hfig,'Cluster',Cluster); end -function UpdateFishDisplay(hfig,new_i_fish) +function UpdateFishDisplay(hfig) % loading steps that are not performed at very first initialization stim = getappdata(hfig,'stim'); -M_fish_set = getappdata(hfig,'M_fish_set'); +fishset = getappdata(hfig,'fishset'); -fishset = M_fish_set(new_i_fish); [~, names] = GetStimRegressor(stim,fishset); + global hstimreg; set(hstimreg,'String',['(choose)',(names)]); @@ -1148,25 +1202,57 @@ function popup_loadfullfishmenu_Callback(hObject,~) set(hfishnum,'Value',new_i_fish); end -function LoadFullFish(hfig,new_i_fish,CONST) +function LoadFullFish(hfig,new_i_fish) data_dir=getappdata(hfig,'data_dir'); -if ~exist('CONST','var'), - disp(['loading fish #' num2str(new_i_fish) '...']); - tic - load(fullfile(data_dir,['CONST_F' num2str(new_i_fish) '.mat']),'CONST'); - toc -end -setappdata(hfig,'isfullfish',1); +disp(['loading fish #' num2str(new_i_fish) '...']); + +%% +tic +fishdir = fullfile(data_dir,['CONST_F' num2str(new_i_fish) '_fast.mat']); +load(fishdir,'const'); +load(fishdir,'dimCR'); +CellResp = zeros(dimCR); +num = 0; +nParts = round(dimCR(1)*dimCR(2)/42000000); +for i = 1:nParts, + load(fishdir,['CRZt_' num2str(i)']); + eval(['len = size(CRZt_' num2str(i) ',1);']); + eval(['CellResp(num+1:num+len,:) = CRZt_' num2str(i) ';']); + eval(['CellResp(num+1:num+len,:) = CRZt_' num2str(i) ';']); + num = num+len; +end +% for i = 1:nParts, +% load(fishdir,['CellResp_' num2str(i)']); +% eval(['len = size(CellResp_' num2str(i) ',1);']); +% eval(['CellResp(num+1:num+len,:) = CellResp_' num2str(i) ';']); +% num = num+len; +% end +toc +setappdata(hfig,'CellResp',CellResp); -% load all fields from CONST, with names preserved -names = fieldnames(CONST); % cell of strings +%% load all fields from CONST, with names preserved +names = fieldnames(const); % cell of strings for i = 1:length(names), - setappdata(hfig,names{i},eval(['CONST.',names{i}])); -end -setappdata(hfig,'numcell',length(CONST.CInfo)); + + + % renaming exception!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if strcmp(names{i},'CRAZ'), + setappdata(hfig,'CellRespAvr',const.CRAZ); + elseif strcmp(names{i},'photostate'), + setappdata(hfig,'stim_full',const.photostate); + + + + else + setappdata(hfig,names{i},eval(['const.',names{i}])); + end +end +setappdata(hfig,'numcell',length(const.CInfo)); + UpdateFishData(hfig,new_i_fish); -UpdateFishDisplay(hfig,new_i_fish); +UpdateFishDisplay(hfig); +setappdata(hfig,'isfullfish',1); end function popup_datamenu_Callback(hObject,~) @@ -1178,29 +1264,49 @@ function popup_datamenu_Callback(hObject,~) end function UpdateDatamenu_Direct(hfig,ID) -CRAZ = getappdata(hfig,'CRAZ'); -CRZt = getappdata(hfig,'CRZt'); -FcAvr = getappdata(hfig,'FcAvr'); +CellResp = getappdata(hfig,'CellResp'); Fc = getappdata(hfig,'Fc'); datanames = getappdata(hfig,'datanames'); tlists = getappdata(hfig,'tlists'); -photostate = getappdata(hfig,'photostate'); +stim_full = getappdata(hfig,'stim_full'); +periods = getappdata(hfig,'periods'); +shift = getappdata(hfig,'shift'); +numcell = getappdata(hfig,'numcell'); +stimtypelist = 1:length(tlists);% getappdata(hfig,'stimtypelist'); + +fishset = getappdata(hfig,'fishset'); + +if fishset == 1, + CellRespAvr = getappdata(hfig,'CellRespAvr'); + FcAvr = getappdata(hfig,'FcAvr'); +else + % find averages of different stimuli, and splice together + % in future should allow user to chose stimtypelist to compile average + % of selected stimuli groups + CRavr = []; + for i = 1:stimtypelist, + M = circshift(CellResp(:,tlists{i}),shift,2); + CRavr = horzcat(CRavr,mean(reshape(M,numcell,periods(i),[]),3)); + end + setappdata(hfig,'CellRespAvr',CRavr); +end + %% SWITCH LEFT/RIGHT DIRECTION AGAIN -FcAvr = vertcat(FcAvr(2,:),FcAvr(1,:),FcAvr(3:end,:)); +% % FcAvr = vertcat(FcAvr(2,:),FcAvr(1,:),FcAvr(3:end,:)); Fc = vertcat(Fc(2,:),Fc(1,:),Fc(3:end,:)); % now 1=left, 2=right, 3=forward, 4 = raw left, 5 = raw right %% if ID == 1, - M_0 = CRAZ; + M_0 = CellRespAvr; fictive = FcAvr; elseif ID == 2, - M_0 = CRZt; + M_0 = CellResp; fictive = Fc; else IX = tlists{ID}; - M_0 = CRZt(:,IX); + M_0 = CellResp(:,IX); temp = Fc(:,IX); % normalize fictive channels, in 2 sets for k = 1:3,%size(F,1), @@ -1211,8 +1317,8 @@ function UpdateDatamenu_Direct(hfig,ID) temp(4:5,:) = (m-min(min(m)))/(max(max(m))-min(min(m))); fictive = temp; end -% stim from photostate -stim = photostate(tlists{ID}); +% stim from stim_full +stim = stim_full(tlists{ID}); setappdata(hfig,'M_0_fluo',M_0); setappdata(hfig,'fictive',fictive); @@ -1411,11 +1517,12 @@ function edit_t_range_Callback(hObject,~) end end -function UpdateTRange(hfig,range) +function UpdateTRange(hfig,range) % this function is incomplete... delete? +% what about stimulus??? + M_0 = GetM_0(hfig); cIX = getappdata(hfig,'cIX'); fictive_0 = getappdata(hfig,'fictive_0'); -% photostate_0 = getappdata(hfig,'photostate_0'); % set M M = M_0(cIX,range); @@ -1423,10 +1530,6 @@ function UpdateTRange(hfig,range) % set fictive fictive = fictive_0(:,range); setappdata(hfig,'fictive',fictive); -% % set photostate -% photostate = photostate_0(:,range); -% setappdata(hfig,'photostate',photostate); - end %% row 2: Operations @@ -1504,17 +1607,17 @@ function popup_ranking_Callback(hObject,~) if length(periods)==1, period = periods{1}; [C,~] = FindCentroid(gIX,M); -else % i_fish>=8, +else %if i_fish==8, tlists = getappdata(hfig,'tlists'); - CRZt = getappdata(hfig,'CRZt'); + CellResp = getappdata(hfig,'CellResp'); IX = tlists{6}; % ptomr_circ - M_ = CRZt(cIX,IX); + M_ = CellResp(cIX,IX); [C,~] = FindCentroid(gIX,M_); periods = getappdata(hfig,'periods'); period = periods{1}+periods{2}; +% else ? + end -% C2 = zscore(C,0,2); -% C_3D = reshape(C2,size(C2,1),period,[]); C_3D_0 = reshape(C,size(C,1),period,[]); C_3D = zscore(C_3D_0,0,2); @@ -1814,9 +1917,57 @@ function PlotRegWithStimMotor(hfig) reg = imNormalize99(regressor); regbar = repmat(reg,[1,1,3]); -subplot(3,1,1);image(stimbar); axis off; title('stimulus'); -subplot(3,1,2);image(regbar); axis off; title('regressor'); -subplot(3,1,3);imagesc(fictive);colormap hot; axis off; title('motor'); +h = subplot(3,1,1);image(stimbar);set(h,'box','on','xtick',[],'ytick',[]); title('stimulus'); +h = subplot(3,1,2);image(regbar); set(h,'box','on','xtick',[],'ytick',[]); title('regressor'); +subplot(3,1,3); + +temp = fictive; +if 1, % plot all 5 lines + fc = vertcat(temp(1,:),temp(3,:),temp(2,:),temp(4,:),temp(5,:)); + + imagesc(fc);colormap gray + set(gca,'YTick',[],'XTick',[]); + set(gcf,'color',[1 1 1]); + set(gca, 'box', 'off') + hold on;axis ij; + + % plot division lines +% for i = 0:3, +% y = i+0.5; +% plot([0.5,length(fictive)+0.5],[y,y],'w','Linewidth',0.5); +% end +% % labels +% names = {'Left','Right','Forward','Raw L','Raw R'}; +% x = -s2*0.05; +% for i = 1:5, +% y = i; +% text(x,y,names{i},'Fontsize',7); +% end + +else % only plot top 3 lines + fc = vertcat(temp(1,:),temp(3,:),temp(2,:)); + imagesc(fc);colormap gray + set(gca,'YTick',[],'XTick',[]); + set(gcf,'color',[1 1 1]); + set(gca, 'box', 'off') + hold on;axis ij; + + % plot division lines + for i = 0:2, + y = i+0.5; + plot([0.5,length(fictive)+0.5],[y,y],'w','Linewidth',0.5); + end + % labels + names = {'Left','Forward','Right'}; + x = -s2*0.07; + for i = 1:3, + y = i; + text(x,y,names{i},'Fontsize',10); + end + +end + +imagesc(fictive);colormap hot; axis off; title('motor'); end function out=imNormalize99(im) @@ -1881,9 +2032,7 @@ function checkbox_flipstim_Callback(hObject,~) stim = getappdata(hfig,'stim'); if regchoice{1}==1, % stim Regressor - i_fish = getappdata(hfig,'i_fish'); - M_fish_set = getappdata(hfig,'M_fish_set'); - fishset = M_fish_set(i_fish); + fishset = getappdata(hfig,'fishset'); [regressors, names] = GetStimRegressor(stim,fishset); regressor = regressors(regchoice{2}).im; @@ -1928,8 +2077,8 @@ function checkbox_flipstim_Callback(hObject,~) fliptrans = GetPhotoTrans(flipstim); % transition e.g.: - % photostate: 2 3 1 3 3 0 0 1 1 0 3 2 0 2 2 1 - % phototrans: 6 11 13 7 15 12 0 1 5 4 3 14 8 2 10 9 + % stimulus : 2 3 1 3 3 0 0 1 1 0 3 2 0 2 2 1 + % transitionID: 6 11 13 7 15 12 0 1 5 4 3 14 8 2 10 9 IX = []; targetstates = unique(fliptrans,'stable'); @@ -2129,11 +2278,11 @@ function pushbutton_IterCentroidRegression_Callback(hObject,~) %% row 3: (TBD) -function checkbox_dataFR_Callback(hObject,~) -hfig = getParentFigure(hObject); -SetDataFR(hfig,get(hObject,'Value')); -RefreshFigure(hfig); -end +% function checkbox_dataFR_Callback(hObject,~) % obsolete !!! +% hfig = getParentFigure(hObject); +% SetDataFR(hfig,get(hObject,'Value')); +% RefreshFigure(hfig); +% end %% ----- tab four ----- (Clustering etc.) @@ -2728,6 +2877,7 @@ function pushbutton_corrplot_Callback(hObject,~) M = getappdata(hfig,'M'); [C,~] = FindCentroid(gIX,M); coeffs = corr(C');%corr(C(1,:)',C(2,:)') +figure('Position',[1000,200,500,500]); CorrPlot(coeffs); end @@ -2744,8 +2894,8 @@ function CorrPlot(coeffs) RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap -figure('Position',[1000,200,500,500]); image(RGB); axis equal; axis tight; +set(gca,'XTick',1:length(im),'YTick',1:length(im)); for i = 1:size(im,2), % horizontal.. because of image axis for j = 1:size(im,1), @@ -3033,7 +3183,7 @@ function SetDataFR(hfig,dataFR) M_0_fluo = getappdata(hfig,'M_0_fluo'); setappdata(hfig,'M',M_0_fluo(cIX,:)); else % Unchecked - M_0_reg = getappdata(hfig,'M_0_reg'); + M_0_reg = getappdata(hfig,'M_0_reg'); % obsolete !!! setappdata(hfig,'M',M_0_reg(cIX,:)); end global hdataFR; @@ -3201,17 +3351,19 @@ function UpdateClusID(hfig,clusID) figure('Position',[100 400 1300 400]); subplot(1,3,1); CORR = corr(C'); - imagesc(CORR);axis equal;axis tight + CorrPlot(CORR); subplot(1,3,2); dendrogram(tree,numU,'orientation','right','reorder',leafOrder); set(gca,'YDir','reverse'); - colormap(jet); + set(gca,'XTick',[]); +% colormap(jet); subplot(1,3,3); C2 = C(leafOrder,:); CORR2 = corr(C2'); - imagesc(CORR2);axis equal;axis tight + CorrPlot(CORR2); +% imagesc(CORR2);axis equal;axis tight end % sort for uniform colorscale temp = zeros(size(gIX)); @@ -3331,7 +3483,7 @@ function UpdateIndices(hfig,cIX,gIX,numK) setappdata(hfig,'cIX',cIX); setappdata(hfig,'gIX',gIX); if exist('numK','var'), - setappdata(hfig,'numK',numK); + setappdata(hfig,'numK',double(numK)); end else errordlg('dataset too large!') @@ -3368,8 +3520,11 @@ function RefreshFigure(hfig) clrmap = getappdata(hfig,'clrmap'); rankscore = getappdata(hfig,'rankscore'); rankID = getappdata(hfig,'rankID'); - iswrite = (rankID>=2); +isPlotLines = 0; %getappdata(hfig,'isPlotLines'); +isPlotFictive = 1; %getappdata(hfig,'isPlotFictive'); + +isPopout = 0; if isempty(cIX), errordlg('empty set! go back!'); @@ -3384,11 +3539,21 @@ function RefreshFigure(hfig) figure(hfig); % left subplot h1 = axes('Position',[0.05, 0.04, 0.55, 0.83]); +% if isCentroid, +% [C,~] = FindCentroid(gIX,M); +% DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); +% else +% DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); +% end if isCentroid, [C,~] = FindCentroid(gIX,M); - DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); + DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,... + iswrite,isPopout,isPlotLines,isPlotFictive); +% DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,iswrite,ispopout); else - DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); + DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,... + iswrite,isPopout,isPlotLines,isPlotFictive); +% DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite,ispopout); end % right subplot diff --git a/GUIpreload_step1_cleanup.m b/GUIpreload_step1_cleanup.m new file mode 100644 index 0000000..f87c089 --- /dev/null +++ b/GUIpreload_step1_cleanup.m @@ -0,0 +1,268 @@ +%%% Loading - Step 1. +% Summary: +% load data from 'cell_resp.stackf' etc, and saved into +% 'FishX_direct_load.mat' (X = number) in the same directory. +% Stimulus parsing is done with function 'StimulusKeyPreprocessing', but +% one is advised to also inspect the results manually. + +% Validation: +% There is an option 'isDiscard50' to discard the noisier 50%, disabled now. +% There is also a section to manually select (draw polygons) regions of +% cells to discard that are anatomical outliers. + +% Format: +% Cell-responses are converted into z-scores and kept that way through out. +% CInfo (cell info for the valid cells) stores the anatomical positions. +% For the rest of parameters saved, see the end of this file. + +% Next: +% For the following Loading steps, only this 'FishX_direct_load.mat' is +% needed, the rest will be cleared from working memory. + +%% +clear all;close all;clc + +%% Set Manaully! + +i_fish = 10; + +M_dir = {'F:\Janelia2014\Fish1_16states_30frames'; + 'F:\Janelia2014\Fish2_20140714_2_4_16states_10frames'; + 'F:\Janelia2014\Fish3_20140715_1_1_16_states_10frames'; + 'F:\Janelia2014\Fish4_20140721_1_8_16states_20frames'; + 'F:\Janelia2014\Fish5_20140722_1_2_16states_30frames'; + 'F:\Janelia2014\Fish6_20140722_1_1_3states_30,40frames'; + 'F:\Janelia2014\Fish7_20140722_2_3_3states_30,50frames'; + 'F:\Janelia2014\Fish8_20141222_2_2_7d_PT_3OMR_shock_lowcut'; + 'F:\Janelia2014\Fish9_20150120_2_1_photo_OMR_prey_blob_blue_cy74_6d_20150120_220356'; + 'F:\Janelia2014\Fish10_20150120_2_2_photo_OMR_prey_blob_blue_cy74_6d_20150120_231917'}; + +fpsec = 1.97; % Hz + +%% load data +datadir = M_dir{i_fish}; +load(fullfile(datadir,'cell_resp_dim.mat')); % 'cell_resp_dim_lowcut.mat' +load(fullfile(datadir,'cell_info.mat')); +cell_resp = read_LSstack_fast_float(fullfile(datadir,'cell_resp.stackf'),cell_resp_dim); +load(fullfile(datadir,'frame_turn.mat')); + +% parse stimulus - should check manually to be sure! +[stimset,stim] = StimulusKeyPreprocessing(frame_turn); + +%% load anatomy +tiffname = fullfile(datadir,'ave.tif'); +info = imfinfo(tiffname,'tiff'); +nPlanes = length(info); +s1 = info(1).Height; +s2 = info(1).Width; +ave_stack = zeros(s1,s2,nPlanes); +for i=1:nPlanes, + ave_stack(:,:,i) = imread(fullfile(datadir,'ave.tif'),i); +end + +% x-y view +im = max(ave_stack,[],3); +out=imNormalize99(im); +anat_yx = repmat(out,[1 1 3]); + +% y-z view +im = squeeze(max(ave_stack,[],2)); +out=imNormalize99(im); +anat_yz = repmat(out,[1 1 3]); + +% x-z view +im = squeeze(max(ave_stack,[],1)); +out=imNormalize99(im); +out = flipud(out'); %%%% empirically necessary... +anat_zx = repmat(out,[1 1 3]); + +dimv_yx = size(anat_yx); +dimv_yz = size(anat_yz); +dimv_zx = size(anat_zx); + +%% compute z-score +cell_resp_z = zscore(cell_resp')'; +nCells = cell_resp_dim(1); + +%% Validating all cells + +isDiscard50 = false; % option to pre-screen +if isDiscard50, + %% Round 1: discard 50% noisy cells based on std of zscore of baseline + + % Compute std of dimest 10% of frames for each cell. + prc = prctile(cell_resp_z,10,2); + STD_full = zeros(nCells,1); + for i = 1:nCells, + ix = find(cell_resp_z(i,:)thr); + [~,I] = sort(STD_full(temp)); + cIX = temp(I); + + % visualize cells to discard + % gIX = (round((1:length(cIX))/1000)+1)'; % option: view sorted index + gIX = round(cIX/1000)+1; % option: view anatomical z index + M = cell_resp_ave_z(cIX,:); + BasicPlotMaps(cIX,gIX,M,cell_info,stim,anat_yx,anat_yz,anat_zx); + + % I_v: index of valid cells + I_v_Holder = ones(1,nCells); + I_v_Holder(cIX) = 0; + I_v = find(I_v_Holder); + + CRAZ = cell_resp_ave_z(I_v,:); + STD = STD_full(I_v); + nCells = length(I_v); + CInfo_0 = cell_info(I_v); + + %% visualize valid cells, sorted by noise + [~,I] = sort(STD); + cIX = I; + gIX = (round((1:length(cIX))/1000)+1)'; % sorted index + M = CRAZ(I,:); + BasicPlotMaps(cIX,gIX,M,CInfo_0,stim,anat_yx,anat_yz,anat_zx); + +else + %% Round 1 alternative: full set + I_v_Holder = ones(1,nCells); + I_v = (1:nCells)'; % valid Index. here not actually screened, I_v is fullset + + cIX = I_v; % 'cell-Index' + gIX = (round((1:length(cIX))/1000)+1)'; % sorted 'group-Index' + M = cell_resp_z(cIX,1:1000); + CInfo_0 = cell_info; + BasicPlotMaps(cIX,gIX,M,CInfo_0,stim,anat_yx,anat_yz,anat_zx); +end + +%% %%%% Round 2: delete anatomically out-of-bound cells by hand +% only execute once (manual choice 1): +cHolder_Anat = []; % collect cIX of all out-of-bound cells + + +% draw cells first +% (choose start and stop index to draw; low numbers ~ ventral) + +%% (optional) step 1: set limit to only plot ventral layer, helps to find outliers within that layer +I_start = 1; +I_stop = round(length(I_v)/10); + +% plot those cells +cIX = I_start:I_stop; +gIX = round(cIX/1000)+1; +figure('Position',[100 0 1300 900]); +numK = round(cIX(end)/1000)+1; +[tot_image, dim_totimage] = DrawClustersOnMap_LSh(CInfo_0,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv','full'); + +% Here manually draw polygons around cells to discard. +% Double click to connect last vertex to first vertex, then double click again within polygon to fix. +% then click again to start drawing the next one +% and when finished, break loop with Ctrl+C... +MaskArray = zeros(dim_totimage(1), dim_totimage(2)); +k_zres = 20; +% draw polygon around outliers on Y-X view +for i = 1:100, % NOMINAL LOOP, break manually (with Ctrl+C, somehow the design didn't work) + h_poly_yx = impoly; + wait(h_poly_yx); % double click to finalize position! + % update finalized polygon in bright color + setColor(h_poly_yx,[0 1 1]); + + IJs = reshape([CInfo_0(I_start:I_stop).center],2,[])'; + A = sub2ind(dimv_yx(1:2),IJs(:,1),IJs(:,2)); + MaskArray = createMask(h_poly_yx); + MaskArray(1:dimv_zx*k_zres+10,:) = []; + B = find(MaskArray); % find indices of pixels within ROI + cIX2 = find(ismember(A,B)); + cHolder_Anat = union(cHolder_Anat,cIX2); + w = waitforbuttonpress; + if w == 1, + break; % this doesn't work ?! + end +end + +%% step 2: plot all cells (including the ones in the ventral layer, don't +I_start = 1; +I_stop = length(I_v); + +% ...here until the end of the cell is the exact dupliate of the last cell... +cIX = I_start:I_stop; +gIX = round(cIX/1000)+1; +figure('Position',[100 0 1300 900]); +numK = round(cIX(end)/1000)+1; +[tot_image, dim_totimage] = DrawClustersOnMap_LSh(CInfo_0,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv','full'); + +% Here manually draw polygons around cells to discard. +% Double click to connect last vertex to first vertex, then double click again within polygon to fix. +% then click again to start drawing the next one +% and when finished, break loop with Ctrl+C... +MaskArray = zeros(dim_totimage(1), dim_totimage(2)); +k_zres = 20; +% draw polygon around outliers on Y-X view +for i = 1:100, % NOMINAL LOOP, break manually (with Ctrl+C, somehow the design didn't work) + h_poly_yx = impoly; + wait(h_poly_yx); % double click to finalize position! + % update finalized polygon in bright color + setColor(h_poly_yx,[0 1 1]); + + IJs = reshape([CInfo_0(I_start:I_stop).center],2,[])'; + A = sub2ind(dimv_yx(1:2),IJs(:,1),IJs(:,2)); + MaskArray = createMask(h_poly_yx); + MaskArray(1:dimv_zx*k_zres+10,:) = []; + B = find(MaskArray); % find indices of pixels within ROI + cIX2 = find(ismember(A,B)); + cHolder_Anat = union(cHolder_Anat,cIX2); + w = waitforbuttonpress; + if w == 1, + break; % this doesn't work ?! + end +end + +%% find extra outliers on top/bottom border of image (can't always get with polygon) +temp = [CInfo_0(:).center]; +XY = reshape(temp',[],nCells)'; +IX = find(XY(:,1)<8 | XY(:,1)> 2040); + +cHolder_Anat = union(cHolder_Anat,IX); + +%% test Plot: all antomy outliers +cIX = cHolder_Anat; +gIX = (1:length(cIX))'; +figure; +DrawClustersOnMap_LSh(CInfo_0,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv','full'); + +%% Clean up anatomy outliers +cIX_Invalid_Anat = I_v(cHolder_Anat); % convert back to index for original full set +I_v_Holder(cIX_Invalid_Anat) = 0; +I_v2 = find(I_v_Holder); + +CR_raw = cell_resp_z(I_v2,:); +nCells = length(I_v2); +CInfo = cell_info(I_v2); + +%% Plot all valid cells to save +cIX = (1:nCells)'; +gIX = round(cIX/1000)+1; +M = CR_raw(cIX,1:1000); +figure; +DrawClustersOnMap_LSh(CInfo,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv','full'); + +%% Save mat files +CInfo_full = cell_info; % save in next round! not saved in current mats!!!!!! + +temp = fullfile(datadir,['Fish' num2str(i_fish) '_full_extrainfo.mat']); +save(temp,'cHolder_Anat','cIX_Invalid_Anat','I_v','I_v2','CInfo_full','-v7.3'); % CInfo_full not saved in mats before 4/24/15!! %% 'STD_full' saved before + +temp = fullfile(datadir,['Fish' num2str(i_fish) '_direct_load.mat']); +varList = {'CR_raw','nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn'}; % used to have 'periods' handcoded +save(temp,varList{:},'-v7.3'); + +%% (if needed) +% save(temp,'somethingelse','-append'); + +%% next run step 2: detrend + + diff --git a/GUIpreload_step2_detrend.m b/GUIpreload_step2_detrend.m new file mode 100644 index 0000000..e17bdc8 --- /dev/null +++ b/GUIpreload_step2_detrend.m @@ -0,0 +1,66 @@ +clear all;close all;clc + +%% manually set directory, number of cores, range of fish (ID) to process + +M_dir = {'F:\Janelia2014\Fish1_16states_30frames'; + 'F:\Janelia2014\Fish2_20140714_2_4_16states_10frames'; + 'F:\Janelia2014\Fish3_20140715_1_1_16_states_10frames'; + 'F:\Janelia2014\Fish4_20140721_1_8_16states_20frames'; + 'F:\Janelia2014\Fish5_20140722_1_2_16states_30frames'; + 'F:\Janelia2014\Fish6_20140722_1_1_3states_30,40frames'; + 'F:\Janelia2014\Fish7_20140722_2_3_3states_30,50frames'; + 'F:\Janelia2014\Fish8_20141222_2_2_7d_PT_3OMR_shock_lowcut'; + 'F:\Janelia2014\Fish9_20150120_2_1_photo_OMR_prey_blob_blue_cy74_6d_20150120_220356'; + 'F:\Janelia2014\Fish10_20150120_2_2_photo_OMR_prey_blob_blue_cy74_6d_20150120_231917'}; + +poolobj=parpool(4); + +range_fish = 10; %1:8 + +%% +for i_fish = range_fish, + disp(['i_fish = ', num2str(i_fish)]); + + datadir = M_dir{i_fish}; + file = fullfile(datadir,['Fish' num2str(i_fish) '_direct_load.mat']); + load(file,'CR_raw'); + % varList = {'CR_raw','nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn','periods'}; + + %% detrend + CR_dtr = zeros(size(CR_raw)); + tmax=size(CR_raw,2); + nCells=size(CR_raw,1); + + tic + parfor i=1:nCells, + cr = CR_raw(i,:); + crd = 0*cr; + for j=1:100:tmax, + if j<=150, + tlim1 = 1; + tlim2 = 300; + elseif j>tmax-150, + tlim1 = tmax-300; + tlim2 = tmax; + else + tlim1 = j-150; + tlim2 = j+150; + end + crr = cr(tlim1:tlim2); + crd(max(1,j-50):min(tmax,j+50)) = prctile(crr,15); + end + if mod(i,100)==0, + disp(num2str(i)); + end + CR_dtr(i,:) = zscore(cr-crd); + end + toc + + CR_dtr = single(CR_dtr); + + %% save into .mat + save(file,'CR_dtr','-append'); +end + +%% +delete(poolobj); diff --git a/GUIpreload_step3.5_PoolClustersFromAllFish.m b/GUIpreload_step3.5_PoolClustersFromAllFish.m new file mode 100644 index 0000000..6ddcbc9 --- /dev/null +++ b/GUIpreload_step3.5_PoolClustersFromAllFish.m @@ -0,0 +1,113 @@ +% Pool from all fish +clear all;close all;clc + +data_dir = 'F:\Janelia2014'; + +global VAR; +load('VAR_current.mat','VAR'); +%% +numfish = 8; +CONSTs = cell(1,numfish); +for i_fish = 1:numfish, + disp(['i_fish = ' num2str(i_fish)]); + + % for new partitioned data, consistant with GUI_050415, but not tested + fishdir = fullfile(data_dir,['CONST_F' num2str(new_i_fish) '_fast.mat']); + load(fishdir,'const'); + load(fishdir,'dimCR'); + CellResp = zeros(dimCR); + num = 0; + nParts = round(dimCR(1)*dimCR(2)/42000000); + for i = 1:nParts, + load(fishdir,['CRZt_' num2str(i)']); + eval(['len = size(CRZt_' num2str(i) ',1);']); + eval(['CellResp(num+1:num+len,:) = CRZt_' num2str(i) ';']); + eval(['CellResp(num+1:num+len,:) = CRZt_' num2str(i) ';']); + num = num+len; + end + names = fieldnames(const); % cell of strings + for i = 1:length(names), + % renaming exception!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +% if strcmp(names{i},'CRAZ'), +% setappdata(hfig,'CellRespAvr',const.CRAZ); +% elseif strcmp(names{i},'photostate'), +% setappdata(hfig,'stim_full',const.photostate); + else + setappdata(hfig,names{i},eval(['const.',names{i}])); + end + end + % add one variable + CONSTs{i_fish}.numcell = length(const.CInfo); + +%% + % load selected clusters + CONSTs{i_fish}.clus = []; + CIX = []; + for i = 1, % just one for now... + if length(VAR(i_fish).ClusGroup{1})-i+1>0, + % copy from VAR + clus = VAR(i_fish).ClusGroup{1}(i); + CONSTs{i_fish}.clus(i).cIX = clus.cIX; + CONSTs{i_fish}.clus(i).gIX = clus.gIX; + CONSTs{i_fish}.clus(i).name = clus.name; + CONSTs{i_fish}.clus(i).numK = clus.numK; + CIX = [CIX; clus.cIX]; + end + end + % copy from CONST + CIX = unique(CIX); + CONSTs{i_fish}.CIX = CIX; + CONSTs{i_fish}.CInfo = const.CInfo(CIX); % CIF + CONSTs{i_fish}.CRAZ = const.CRAZ(CIX,:); + CONSTs{i_fish}.CRZt = CellResp(CIX,:); + + %% OLD FORMAT without partitioning +% % load general info +% load(['CONST_F' num2str(i_fish) '.mat'],'CONST'); +% % CONST: {'ave_stack','anat_yx','anat_yz','CInfo','periods','shift',... +% % 'CRAZ','CRZt','dshift_fc','FcAvr','Fc','photostate','tlists','datanames'}; +% %% +% % load fields from CONST, with names preserved +% names = {'anat_yx','anat_yz','anat_zx','periods','shift','dshift','FcAvr','Fc','photostate','tlists','datanames'}; +% +% for i = 1:length(names), +% eval(['CONSTs{i_fish}.',names{i},' = CONST.',names{i},';']); +% end +% +% % add one variable +% CONSTs{i_fish}.numcell = length(CONST.CInfo); +% %% +% % load selected clusters +% CONSTs{i_fish}.clus = []; +% CIX = []; +% for i = 1, % just one for now... +% if length(VAR(i_fish).ClusGroup{1})-i+1>0, +% % copy from VAR +% clus = VAR(i_fish).ClusGroup{1}(i); +% CONSTs{i_fish}.clus(i).cIX = clus.cIX; +% CONSTs{i_fish}.clus(i).gIX = clus.gIX; +% CONSTs{i_fish}.clus(i).name = clus.name; +% CONSTs{i_fish}.clus(i).numK = clus.numK; +% CIX = [CIX; clus.cIX]; +% end +% end +% % copy from CONST +% CIX = unique(CIX); +% CONSTs{i_fish}.CIX = CIX; +% CONSTs{i_fish}.CInfo = CONST.CInfo(CIX); % CIF +% CONSTs{i_fish}.CRAZ = CONST.CRAZ(CIX,:); +% CONSTs{i_fish}.CRZt = CONST.CRZt(CIX,:); +end +%% +clear CONST; +temp = whos('CONSTs'); +if [temp.bytes]>2*10^9, + save('CONSTs_current.mat','CONSTs','-v7.3'); +else + save('CONSTs_current.mat','CONSTs'); +end + +%% optional +timestamp = datestr(now,'mmddyy_HHMM'); +currentfolder = pwd; +save([currentfolder '\arc mat\CONSTs_' timestamp '.mat'],'CONSTs'); diff --git a/GUIpreload_step3_align.m b/GUIpreload_step3_align.m new file mode 100644 index 0000000..864e55d --- /dev/null +++ b/GUIpreload_step3_align.m @@ -0,0 +1,298 @@ +%%%%%%%%%%%%%%% +%% +% clear all;close all; +% varList = {'CR_dtr','nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn','periods'}; + +% % or +varList = {'nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn'}; + +%% +M_dir = {'F:\Janelia2014\Fish1_16states_30frames'; + 'F:\Janelia2014\Fish2_20140714_2_4_16states_10frames'; + 'F:\Janelia2014\Fish3_20140715_1_1_16_states_10frames'; + 'F:\Janelia2014\Fish4_20140721_1_8_16states_20frames'; + 'F:\Janelia2014\Fish5_20140722_1_2_16states_30frames'; + 'F:\Janelia2014\Fish6_20140722_1_1_3states_30,40frames'; + 'F:\Janelia2014\Fish7_20140722_2_3_3states_30,50frames'; + 'F:\Janelia2014\Fish8_20141222_2_2_7d_PT_3OMR_shock_lowcut'; + 'F:\Janelia2014\Fish9_20150120_2_1_photo_OMR_prey_blob_blue_cy74_6d_20150120_220356'; + 'F:\Janelia2014\Fish10_20150120_2_2_photo_OMR_prey_blob_blue_cy74_6d_20150120_231917'}; + +M_stimset = [1, 1, 1, 1, 1, 1, 1, 2, 2, 2]; + +save_dir = 'F:\Janelia2014'; + +%% MANUAL +for i_fish = 10, %:8, + disp(num2str(i_fish)); + %% load data + datadir = M_dir{i_fish}; + load(fullfile(datadir,['Fish' num2str(i_fish) '_direct_load.mat']),varList{:}); + + %% ONLY RUN ONCE!!!!!!!!! + % add 1 frame at start: + CR_dtr = horzcat(CR_dtr(:,1),CR_dtr); + + frame_turn = vertcat(frame_turn(1,:),frame_turn); + + % correction of an error in Fish #1 + if i_fish == 1, + CR_dtr = horzcat(CR_dtr(:,1:20),CR_dtr); + end + + %% index processing + if M_stimset(i_fish)==1, % fish 1-7 + M_period = {480,160,160,320,480,280,300}; %,{120,150,360}}; + periods = M_period{i_fish}; + period = periods; + + nrep = floor(size(CR_dtr,2)/period)-1; + dshift = 2; + shift = period+dshift; + + IX_all = 1+shift:period*nrep+shift; + CellResp = CR_dtr(:,IX_all); % old name: CRZt, 'Cell Responses Zscore trimmed' + + IX_avr = 1+shift:period+shift; + nCells = size(CR_dtr,1); + CRAZ = mean(reshape(CR_dtr(:,IX_all),nCells,period,[]),3); % 'Cell Response Average Zscore' + + datanames = {'rep average','all reps','rep #1','rep #2','last rep'}; + ix_avr = 1:period; + ix_all = 1:period*nrep; + tlists = {ix_avr, ix_all}; + for i = [1,2,nrep], + ix = 1+period*(i-1):period*i; + tlists = [tlists, ix]; + end + + stim_full = frame_turn(17,IX_all); + stim_full = round(stim_full); + + else % multiple protocols + % process raw stimulus code + [stimset,stim_full] = StimulusKeyPreprocessing(frame_turn); + + % find time-lists = list of time-frames numbers for a certain stimulus + nTypes = length(stimset); + tlists_raw = cell(1,nTypes+1); % time-lists, i.e. frame indices + IX_all = []; + for i_type = 1:nTypes, + nSets = length(stimset(i_type).starts); + IX = []; + for i_set = 1:nSets, + IX = horzcat(IX,stimset(i_type).starts(i_set):stimset(i_type).stops(i_set)); + end + tlists_raw{i_type} = IX; + IX_all = horzcat(IX_all,IX); + end + tlists_raw{nTypes+1} = IX_all; + + % This is the main data to store and load to GUI + CellResp = CR_dtr(:,IX_all); % old name: CRZt, 'Cell Responses Zscore trimmed' + + %% get tlists corresponding to IX_all (tlists_raw ~ 1:totalframe#) + tlists = cell(1,nTypes+1); % time-lists, i.e. frame indices + for i_type = 1:nTypes, + [~,tlists{i_type}] = intersect(IX_all,tlists_raw{i_type}); + end + + %% find average + CRAZ = []; % 'Cell Response Average Zscore' + stimtypelist = 1:nTypes; + numcell = size(CR_dtr,1); + for i = stimtypelist, + M = circshift(CR_dtr(:,tlists_raw{i}),shift,2); + CRAZ = horzcat(CRAZ,mean(reshape(M,numcell,periods(i),[]),3)); + end + + %% variables to save later in struct + periods = [stimset.period]; + datanames = {'rep average','all reps',stimset.name}; + + dshift = 2; + shift = -dshift; % circshift fluo left by 2 frames + + end + + %% prepare fictive data + rows = [7,8,9,13,14]; + F = frame_turn(:,rows)'; + + F(2,:) = -F(2,:); + Fc = F(:,IX_all); + + % find averages + if M_stimset(i_fish)==1, % fish 1-7 + FcAvr = mean(reshape(Fc,length(rows),period,[]),3); + else + FcAvr = []; % 'Cell Response Average Zscore' + stimtypelist = 1:nTypes; + numcell = size(CR_dtr,1); + for i = stimtypelist, + avr = mean(reshape(F(:,tlists_raw{i}),length(rows),periods(i),[]),3); + FcAvr = horzcat(FcAvr,avr); + end + end + + % normalizations + for i = 1:3, + m = FcAvr(i,:); + FcAvr(i,:) = (m-min(m))/(max(m)-min(m)); + m = Fc(i,:); + Fc(i,:) = (m-min(m))/(max(m)-min(m)); + end + m = FcAvr(4:5,:); + FcAvr(4:5,:) = (m-min(min(m)))/(max(max(m))-min(min(m))); + m = Fc(4:5,:); + Fc(4:5,:) = (m-min(min(m)))/(max(max(m))-min(min(m))); + + + %% compile CONST + CONST = []; + names = {'ave_stack','anat_yx','anat_yz','anat_zx','CInfo','periods','shift','CellResp','CRAZ','dshift','Fc','FcAvr','stim_full','stimset','tlists','datanames'}; + for i = 1:length(names), % use loop to save variables into fields of CONST + eval(['CONST.',names{i},'=', names{i},';']); + end + + %% and save + %%% old method: + % temp = whos('CONST'); + % if [temp.bytes]<2*10^9, + % save(['CONST_F' num2str(i_fish) '.mat'],'CONST'); + % beep; + % else + % save(['CONST_F' num2str(i_fish) '.mat'],'CONST','-v7.3'); + % beep; + % end + + %%% new method with partitioning of main data + newfishdir = fullfile(save_dir,['CONST_F' num2str(i_fish) '_fast.mat']); + const = CONST; + const = rmfield(const,'CellResp'); + dimCR = size(CONST.CellResp); + save(newfishdir,'const','-v6'); + save(newfishdir,'dimCR','-v6','-append'); + + nParts = round(numel(CONST.CellResp)/42000000); + t = zeros(1,nParts); + ix = round(linspace(1,size(CONST.CellResp,1),nParts+1)); + for i = 1:nParts, + disp(num2str(i)); + eval(['CellResp_' num2str(i) '= CONST.CellResp(ix(i):ix(i+1),:);']); + save(newfishdir,['CellResp_' num2str(i)],'-v6','-append'); + end + + + %% OR + +end + + %% initialize VAR (once) + + +%% Initialize VARS % outdated? Clusgroup? + +nCells = length(CONST.CInfo); + + cIX = (1:nCells)'; + i = 1; + VAR(i_fish).Class(i).name = 'all processed'; + VAR(i_fish).Class(i).cIX = cIX; + VAR(i_fish).Class(i).gIX = ones(length(cIX),1); + VAR(i_fish).Class(i).numel = nCells; + VAR(i_fish).Class(i).numK = 1; + VAR(i_fish).Class(i).datatype = 'std'; + + cIX = (1:100:nCells)'; + VAR(i_fish).ClusGroup{1,1}.name = 'test'; + VAR(i_fish).ClusGroup{1,1}.cIX = cIX; + VAR(i_fish).ClusGroup{1,1}.gIX = ones(length(cIX),1); + VAR(i_fish).ClusGroup{1,1}.numel = length(cIX); + VAR(i_fish).ClusGroup{1,1}.numK = 1; + + %% + cIX = (1:10:nCells)'; + VAR(i_fish).ClusGroup{1,1}(12).name = '1/10 of all'; + VAR(i_fish).ClusGroup{1,1}(12).cIX = cIX; + VAR(i_fish).ClusGroup{1,1}(12).gIX = ones(length(cIX),1); + VAR(i_fish).ClusGroup{1,1}(12).numel = length(cIX); + VAR(i_fish).ClusGroup{1,1}(12).numK = 1; + + %% varience/std for reps for each cell +% if i_fish==2 || i_fish==3 || i_fish==6, +% period_real = period/2; +% else +% period_real = period; +% end +% nrep_real = floor((size(CR,2)-shift)/period_real); +% while period_real*nrep_real+shift>size(CR,2), +% nrep_real = nrep_real-1; +% end +% CRZ_3D = reshape(CRZ(:,1+shift:period_real*nrep_real+shift),nCells,period_real,[]); +% %% updated method, weighing both std between each rep and (summed with) std of 1st half & 2nd half of experiment - 1/8/15 +% % CRZ = CONST.M_array.CellResp; +% % if i_fish==2 || i_fish==3 || i_fish==6, +% % period_real = CONST.M_array.period/2; +% % else +% % period_real = CONST.M_array.period; +% % end +% % CRZ_3D = reshape(CRZ,size(CRZ,1),period_real,[]); +% % divide = round(size(CRZ_3D,3)/2); +% % CRZ_std1 = std(CRZ_3D(:,:,1:divide),0,3); +% % CRZ_std2 = std(CRZ_3D(:,:,divide+1:end),0,3); +% % temp1 = mean(CRZ_std1,2); +% % temp2 = mean(CRZ_std2,2); +% % +% % temp12 = horzcat(temp1,temp2); +% % temp = mean(temp12,2)+std(temp12,0,2); +% % [~,I] = sort(temp); +% % M = temp(I); +% % figure;plot(M) +% % +% % figure;imagesc(CRZ(I,:)) +% % +% % nCells = size(CRZ,1); +% +% %% find low variance / stimulus-locked cells +% CRZ_std = std(CRZ_3D,0,3); +% temp = mean(CRZ_std,2); +% +% % find mean-std thres: 0.5 +% [~,I] = sort(temp); +% M = temp(I); +% figure;plot(M) +% %% +% i_last = length(VAR(i_fish).Class); +% M_perc = [0.025,0.1,0.3]; +% for j = 1:length(M_perc); +% thres = M(round(nCells*M_perc(j))); +% cIX = find(tempthr) = thr+1; % assign cap value to be thr+1 + +xv = 1:thr+1; % x-vector for histogram bins + +% Manual inspection +% figure; +% hist(st,xv); + +% (continue by default) find bins with significantly high counts, i.e. +% whole range/set of raw stimulus keys, store in 'M_range_raw' +% use to manually define standardized 'M_range' +counts = hist(st,xv); +thr_counts = round(length(st)/500); +M_range_raw = xv(counts>thr_counts); +disp(['M_range_raw = ' num2str(M_range_raw)]); % use this to adjust the substitution array, M_range, manually after this inspection!%%%%%%%%%%%%%%%%%%%%%% + +%% Standardized Stimulus Key ----------------------------------------------------------------------------- MANUAL INPUT! +% Define normalized code, e.g. same stimulus under different raw stim keys +% can be unified +M_range = [3,1,3,2, 9,10,11,12, 0,13, 3,14,15, 4,16]; % corresponding to +% raw = [2,3,4,5,12,13,14,15,22,23,30,31,32,99,100] % for Fish 10 + +% 0 = all black; 1 = black/white; 2 = white/black; 3 = all white; 4 = all gray; +% 5 = gray/black; 6 = white/gray; 7 = black/gray; 8 = gray/white. +% OMR baseline = not moving?? grey?? +% 10 = forward grating (very slow, more for calibration) +% 11 = rightward grating +% 12 = leftward grating + + +% half-fields can be displayed directly. other L/R directional cues can be +% displayed by gradient. If not-directional or unspecified, will assign a +% color-code randomly. + +%% correct transient frames +% short method: st_rounded = interp1(M_range,M_range,st,'nearest'); +% but this doesn't work if a frame is different from both the one +% before and the one after. Use loop below to force a manual 'round' +% to the closer one of the 2 neighbors. +temp = diff([0,st]); +switches = find(temp); +trans = switches(diff(switches)==1); + +for i = 1:length(trans), + j = trans(i); + if j>1, + st(j) = st(j-1) + round((st(j)-st(j-1))/abs(st(j+1)-st(j-1))); % round to nearest of st(j-1) and st(j+1) + end +end + +% (loop until all transient ones are corrected) +temp = diff([0,st]); +switches = find(temp); +trans = switches(diff(switches)==1); +count = 0; +while ~isempty(trans), + for i = 1:length(trans), + j = trans(i); + st(j) = st(j-1) + round((st(j)-st(j-1))/abs(st(j+1)-st(j-1))); % round to nearest of st(j-1) and st(j+1) + end + temp = diff([0,st]); + switches = find(temp); + trans = switches(diff(switches)==1); + count = count+1; + if count>5, + disp('count>5??'); + break; + end +end + +%% Standardize stimulus keys for both 'stim' and 'block', based on 'M_range_raw'->'M_range' +% (tricky for cell arrays) custom method: +% make a replacement matrix 'M_replace' to substitude raw stimulus key in 'block_raw' +% to normalized key, stored in new cell array 'block' +M_replace = zeros(1,max(M_range_raw)); +for i = 1:max(M_range_raw), + ix = find(i==M_range_raw); + if ~isempty(ix), + M_replace(i) = M_range(ix); + end +end + +stim = zeros(size(st)); +block = cell(nBlocks,1); +for i = 1:length(M_range_raw), + stim(st==M_range_raw(i)) = M_range(i); + for j = 1:nBlocks, + block{j} = cellfun(@(x) M_replace(x),block_raw{j},'UniformOutput',0); + end +end + +%% plot stimulus-key! + + +%% Stimulus set segmentation, based on 'block' (manual input of protocol sequence) +% reduce 'stim' to a sequence of keys without consecutive repeats +ix_singles = find(diff([0,stim])); % index array to map back to raw indices (frame-number) +singles = stim(ix_singles); % sequence of keys without consecutive repeats. +% '_sg' below stands for 'singles'. Manipulations below in both raw and singles, in parallel. + +jump = 1; % searching for 'jumps' in 'singles' +% Based on the known protocol sequence for each set, look for first +% mismatch in 'singles', i.e. position where the next set begins, stored +% as 'block_change(_sg)' and 'set_start(_sg)'. +block_change = cell(nBlocks,1); +block_change_sg = cell(nBlocks,1); +for I = 1:nBlocks, + nSets = length(block{I}); + block_change{I} = zeros(1,nSets); + block_change_sg{I} = zeros(1,nSets); + for J = 1:nSets, + jump_last = jump; + pattern = block{I}{J}; % protocol sequence for this set + + % this is sort of awkward. Instead of padding the 'already-found' sets with zeros, + % one could also search the cropped array and adjusted the index. + % Now 'comparison' is aligned with stim, works. + textile = repmat(pattern,1,ceil(length(singles)/length(pattern))); % tile, for use below + + template = zeros(size(singles)); + ixs_fill = jump_last:length(singles); + template(ixs_fill) = textile(1:length(ixs_fill)); % crop the 'textile' to fill the space after jump_last + + singles_temp = singles; + if jump_last>1, + singles_temp(1:jump_last-1) = 0; + end + comparison = (singles_temp-template); + % find the position where the next set begins + jump = find((comparison),1,'first'); + + %% contingency plan: if start of the new set is different from + % expected protocol ('block') + count_stuck = 0; + count_stuck2 = 0; + while jump-jump_lastlength(pattern), + count_stuck = 0; + jump_last = jump_last+1; + count_stuck2 = count_stuck2+1; + disp(['count_stuck2=' num2str(count_stuck2)]); + if count_stuck2>10, % arb. thres + disp('count_stuck2>10'); + break; + end + end + + pattern_ = circshift(pattern,-count_stuck,2); + textile = repmat(pattern_,1,ceil(length(singles)/length(pattern_))); + + template = zeros(size(singles)); + ixs_fill = jump_last:length(singles); + template(ixs_fill) = textile(1:length(ixs_fill)); + + singles_temp = singles; + if jump_last>1, + singles_temp(1:jump_last-1) = 0; + end + comparison = (singles_temp-template); + + jump = find((comparison),1,'first'); + + end + %% save + if isempty(jump), % reached end of experiment + %% + % NOTE END OF + % EXPERIMENT???????????????????????????????????????????????????????????????????????? + break; + end + block_change_sg{I}(J) = jump; + block_change{I}(J) = ix_singles(jump); + end +end + +% eliminate empty sets +for I = 1:nBlocks, + nSets = length(block_change{I}); + for J = 1:nSets, + if block_change{I}(J) == 0, + block_change{I}(J:end) = []; + block_change_sg{I}(J:end) = []; + if length(block{I})>= J, + block{I}(J:end) = []; + end + break; + end + end +end + +% block_change is basically start of next set. shift by one set to get 'set_start'. +set_start = cell(nBlocks,1); +set_start_sg = cell(nBlocks,1); +set_stop = cell(nBlocks,1); +set_stop_sg = cell(nBlocks,1); +for I = 1:nBlocks, + nSets = length(block{I}); + set_start{I} = zeros(1,nSets); + set_start_sg{I} = zeros(1,nSets); + if I == 1, + set_start{I}(1) = 1; + set_start_sg{I}(1) = 1; + else + set_start{I}(1) = block_change{I-1}(end); + set_start_sg{I}(1) = block_change_sg{I-1}(end); + end + set_start{I}(2:end) = block_change{I}(1:end-1); + set_start_sg{I}(2:end) = block_change_sg{I}(1:end-1); + set_stop{I} = block_change{I}-1; + set_stop_sg{I} = block_change_sg{I}-1; +end + +%% Visualize current segmentation of blocks +if isPlotting, + figure; + hold on; + plot(stim) + for I = 1:nBlocks, + for J = 1:length(block_change{I}), + x = block_change{I}(J); + plot([x,x],[0,20],'r--') + end + end +end +%% Find shift-corrections for each stim-set, used for averaging later + + +%% Find periods for each set-stype, and organize info by type into 'stimset'. +for i_ss = 1:length(stimset), + % find empty sets (for this particular experiment) and delete ij's in 'stimset' + for i_SetRep = 1:size(stimset(i_ss).ij,1), + I = stimset(i_ss).ij(i_SetRep,1); + J = stimset(i_ss).ij(i_SetRep,2); + if I>length(block) || J>length(block{I}), + stimset(i_ss).ij(i_SetRep:end,:) = []; + break; + end + end + nSetReps = length(stimset(i_ss).ij); + + %% get period from first set of block + I = stimset(i_ss).ij(1,1); + J = stimset(i_ss).ij(1,2); + pattern = block{I}{J}; % protocol sequence (singles) in normalized code + stimset(i_ss).pattern = pattern; + + % find periodicity + if length(pattern)==1, % exception, e.g. 'Spontaneous' set, single stimlus key held for full duration of set + if J < length(block{I}), + next_start = set_start{I}(J+1); + elseif J == length(block{I}), + next_start = set_start{I+1}(1); + else + next_start = length(stim); + end + stimset(i_ss).period = next_start - set_start{I}(J); + else % find periodicity based on stimlus key repetition + for i = 1:length(pattern), + if length(find(pattern==pattern(i)))==1, % use a stimulus that appears only once per period + % find period + start_s = set_start_sg{I}(J) + i -1; + stop_s = start_s + find(singles(1,start_s+1:end)==pattern(i),1,'first'); + stimset(i_ss).period = ix_singles(stop_s)-ix_singles(start_s); + break; + end + end + end + + %% + % stimset(i_ss).period % assigned above + + stimset(i_ss).rawstarts = zeros(1,nSetReps); + stimset(i_ss).rawstops = zeros(1,nSetReps); + stimset(i_ss).starts = zeros(1,nSetReps); + stimset(i_ss).stops = zeros(1,nSetReps); + stimset(i_ss).nReps = zeros(1,nSetReps); + for i_SetRep = 1:nSetReps, + I = stimset(i_ss).ij(i_SetRep,1); + J = stimset(i_ss).ij(i_SetRep,2); + stimset(i_ss).rawstarts(i_SetRep) = set_start{I}(J); + stimset(i_ss).rawstops(i_SetRep) = block_change{I}(J)-1; + + if length(pattern)==1, % exception, e.g. 'Spontaneous' set, single stimlus key held for full duration of set + stimset(i_ss).starts(i_SetRep) = stimset(i_ss).rawstarts(i_SetRep); + stimset(i_ss).stops(i_SetRep) = stimset(i_ss).rawstops(i_SetRep); + stimset(i_ss).nReps(i_SetRep) = 1; + else + % circshift to find intact periods + pattern = stimset(i_ss).pattern; + + ix_start = set_start_sg{I}(J); + ix_stop = set_stop_sg{I}(J); + + % find 'shift' (shift = 1 means starting from 2nd in 'singles') + segment = singles(ix_start:ix_start+length(pattern)-1); + shift = 0; + while ~isequal(segment,pattern), + shift = shift+1; + segment = singles(ix_start+shift:ix_start+length(pattern)-1+shift); + end + stimset(i_ss).starts(i_SetRep) = ix_singles(ix_start+shift); + % round to period within this set + nReps = floor((ix_stop-ix_start-shift+1)/length(pattern)); + stimset(i_ss).nReps(i_SetRep) = nReps; + + ix = stimset(i_ss).starts(i_SetRep) + nReps*stimset(i_ss).period - 1; + if ix>stimset(i_ss).rawstops(i_SetRep), + nReps = nReps-1; + end + stimset(i_ss).stops(i_SetRep) = stimset(i_ss).starts(i_SetRep) + nReps*stimset(i_ss).period - 1; + end + end +end + +% discard very first period because of potential artifacts at beginning of +% entire experiment +if stimset(1).nReps(1)>1, + stimset(1).nReps(1) = stimset(1).nReps(1) - 1; + stimset(1).starts(1) = stimset(1).starts(1) + stimset(1).period; +end + + +end \ No newline at end of file diff --git a/arrows.jpg b/arrows.jpg new file mode 100644 index 0000000..6d1147b Binary files /dev/null and b/arrows.jpg differ diff --git a/arrows_wide.jpg b/arrows_wide.jpg new file mode 100644 index 0000000..9438967 Binary files /dev/null and b/arrows_wide.jpg differ diff --git a/first version/AccessLocalFunctions.m b/first version/AccessLocalFunctions.m new file mode 100644 index 0000000..1953044 --- /dev/null +++ b/first version/AccessLocalFunctions.m @@ -0,0 +1,12 @@ + + +f = []; +for i = 1:length(fcns), + eval(['f.' char(fcns{i}) ' = fcns{i};']); +end + +%% test function 'FindCentroid' +M = getappdata(hfig,'M'); +gIX = getappdata(hfig,'gIX'); +C_test = f.FindCentroid(gIX,M); +figure;imagesc(C_test) \ No newline at end of file diff --git a/first version/BasicPlotMaps.m b/first version/BasicPlotMaps.m new file mode 100644 index 0000000..9196801 --- /dev/null +++ b/first version/BasicPlotMaps.m @@ -0,0 +1,261 @@ +function BasicPlotMaps(cIX,gIX,M,CInfo,photostate,anat_yx,anat_yz,anat_zx) %(cIX,gIX,cIX_0,M_0,CIF,numK,M) +% M = M_0(cIX_0(cIX),:); +numK = length(unique(gIX)); +%% +dataFR = 1; +figure('Position',[100 50 1350 900]);%,'DeleteFcn',@closefigure_Callback); +h1 = axes('Position',[0.05, 0.06, 0.5, 0.85]); % left ~subplot +BasicDrawClusters(h1,M,gIX,dataFR,numK,photostate); + +h2 = axes('Position',[0.58, 0.06, 0.40, 0.85]); % right ~subplot +DrawClustersOnMap_LSh(CInfo,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv'); +% set(gcf,'PaperUnits', 'inches', 'PaperSize', [12, 6]) + +end + +function BasicDrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite) +pos = get(gca,'Position'); +barratio = 0.02; +numK = double(max(double(numK),double(max(gIX)))); + +% sort traces by index +im = M; +[nLines,nFrames] = size(im); +[~,I] = sort(gIX); +im = im(I,:); + +%% convert imagesc effect into rgb matrix 'RGB' +if dataFR, % is drawing cell-response fluorescence + cmap = gray(64); + im = AutoScaleImage0to1(im); % scaling to min/max 0/1 + minlim = 0; + maxlim = 1; +else % create blue-white-red map, for regression results + cmap = zeros(64,3); + cmap(:,1) = [linspace(0,1,32), linspace(1,1,32)]; + cmap(:,2) = [linspace(0,1,32), linspace(1,0,32)]; + cmap(:,3) = [linspace(1,1,32), linspace(1,0,32)]; + minlim = -1; %min(min(im)); + maxlim = 1; %max(max(im)); +end +RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap + +%% add vertical color code +if exist('clrmap','var'), + if strcmp(clrmap,'jet'), + temp = flipud(jet(numK)); + else % 'hsv' + temp = hsv(round(numK*1.1)); + end +else % 'hsv' + temp = hsv(round(numK*1.1)); +end +cmap2 = vertcat(temp(1:numK,:),[0,0,0]); % extend colormap to include black +bwidth = max(round(nFrames/30),1); + +idx = gIX(I); +ix_div = [find(diff(idx));length(idx)]; + +bars = ones(nLines,bwidth); +if numK>1, + bars(1:ix_div(1),:) = idx(ix_div(1)); + for i = 2:length(ix_div), + % paint color + bars(ix_div(i-1)+1:ix_div(i),:) = idx(ix_div(i)); + end +end +im_bars = reshape(cmap2(bars,:),[size(bars),3]); +% add a white division margin +div = ones(nLines,round(bwidth/2),3); +% put 3 parts together +im = horzcat(RGB,div,im_bars); + +%% plot figure +% h_fig = figure('Position',[100 0 1300 900],'Name',['round ' num2str(i_step)]); %[left, bottom, width, height] +[s1,s2,s3] = size(im); + +hold off; +set(h1,'Position',[pos(1),pos(2)+pos(4)*barratio*3,pos(3),pos(4)*(1-4*barratio)]); +if s1<30, + temp = im; + im = ones(30,s2,s3); + im(1:s1,:,:) = temp; +end +image(im); +set(gca, 'box', 'off') + +hold on; + +% plot cluster division lines +if numK>1, + for i = 1:length(ix_div),% = numK-1, + y = ix_div(i)+0.5; + plot([0.5,s2+0.5],[y,y],'k','Linewidth',0.5); +% plot([1,s2*kx],[y,y],'k','Linewidth',0.5); + end +end + +ylabel(['ROI number: ' num2str(s1)]); +set(gca,'YTick',[],'XTick',[]); +set(gcf,'color',[1 1 1]); + +colormap gray; + +% write in number label +x = s2*1.003; +y_last = 0; +for i = 1:length(ix_div), + % avoid label crowding + margin = 0.015*s1; + y0 = ix_div(i)+0.5; + y = max(y_last+margin,y0); + if ilow_high(2)) = low_high(2); +im = mat2gray(temp); +end + +function RGB = ImageToRGB(im,cmap,minlim,maxlim) +L = size(cmap,1); +ix = round(interp1(linspace(minlim,maxlim,L),1:L,im,'linear','extrap')); +RGB = reshape(cmap(ix,:),[size(ix) 3]); % Make RGB image from scaled. +end + +function DrawBars(pos,barratio,nFrames,barlength,photostate) % horizontal stimulus bar +% pos = get(gca,'Position'); % [0.05, 0.06, 0.5, 0.85] + +% global htop hbottom; +barheight = 100; +stimheight = 2*barheight; + +m1 = 0.3*ones(barheight,length(photostate)); +m2 = m1; % bottom half +x = photostate; + +% 2,3,4,5,12,13,14,15,99 +% 3,1,3,2,4 ,10,11,12, + +if max(x)<=3, + % 0 = all black; 1 = black/white; 2 = white/black; 3 = all white; (4 = all gray;) + % 5 = gray/black; 6 = white/gray; 7 = black/gray; 8 = gray/white. + + % top (~ projection right) + m1(:,x==0 | x==2 | x==5) = 0; % black + m1(:,x==4 | x==6 | x==7) = 0.5; % grey + m1(:,x==1 | x==3 | x==8) = 1; % white + % bottom (~ projection left) + m2(:,x==0 | x==1 | x==7) = 0; + m2(:,x==4 | x==5 | x==8) = 0.5; + m2(:,x==2 | x==3 | x==6) = 1; + + temp = repmat(vertcat(m1,m2),[1 1 3]); % 3 color layers in dim3 +else + %% + % 2 = white/white + % 3 = black/white + % 4 = white/white + % 5 = white/black + % + % 12 = gray + % 13 = forward grating (very slow, more for calibration) + % 14 = rightward grating + % 15 = leftward grating (* I hope I got right/left correct - if you think it's flipped it probably is) + % + % 99 = spontaneous + % top (~ projection right) + + m1(:,x==3) = 0; % black + m1(:,x==2 | x==4 | x==5) = 1; % white + % bottom (~ projection left) + m2(:,x==5) = 0; + m2(:,x==2 | x==3 | x==4) = 1; + + temp = repmat(vertcat(m1,m2),[1 1 3]); % 3 color layers in dim3 + %% + temp(:,x==14,1) = 1; % red + temp(:,x==13,2) = 1; % green + temp(:,x==15,3) = 1; % blue + + % grey + temp(:,x==12,1) = 0.8; + temp(:,x==12,2) = 0.8; + temp(:,x==12,3) = 0.8; + + % yellow + temp(:,x==99,1) = 1; + temp(:,x==99,2) = 1; + temp(:,x==99,3) = 0.6; +end + +n = ceil(nFrames/length(photostate)); % tile, to size or bigger +temp = repmat(temp,1,n); +stimbar = ones(stimheight,barlength,3); +stimbar(:,1:nFrames,:) = temp(:,1:nFrames,:); % crop to size + +% if isempty(htop), +axes('Position',[pos(1),pos(2)+pos(4)*(1-barratio),pos(3),pos(4)*barratio]); +% elseif ~ishandle(htop), +% htop = axes('Position',[pos(1),pos(2)+pos(4)*(1-barratio),pos(3),pos(4)*barratio]); +% else axes(htop) +% end +image(stimbar);axis off;hold on; +% plot top border +plot([1,length(photostate)],[1,1],'k','Linewidth',0.5); + +% if isempty(hbottom), +axes('Position',[pos(1),pos(2)+pos(4)*barratio*2,pos(3),pos(4)*barratio]); +% elseif ~ishandle(hbottom), +% hbottom = axes('Position',[pos(1),pos(2),pos(3),pos(4)*barratio]); +% else axes(hbottom) +% end + +image(stimbar);axis off;hold on; +% plot bottom border +plot([1,length(photostate)],[stimheight,stimheight],'k','Linewidth',0.5); +%% +end diff --git a/Conv_test.m b/first version/Conv_test.m similarity index 100% rename from Conv_test.m rename to first version/Conv_test.m diff --git a/DA_2_Solution.m b/first version/DA_2_Solution.m similarity index 100% rename from DA_2_Solution.m rename to first version/DA_2_Solution.m diff --git a/first version/DrawClusters.m b/first version/DrawClusters.m new file mode 100644 index 0000000..964da88 --- /dev/null +++ b/first version/DrawClusters.m @@ -0,0 +1,210 @@ +function DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite) +pos = get(gca,'Position'); +barratio = 0.03; + +%% Prepare cluster data +% down-sample +displaymax = 1000; +numcell = size(M,1); +if numcell > displaymax, + skip = round(numcell/displaymax); + M = M(1:skip:end,:); + gIX = gIX(1:skip:end,:); +end + +numK = double(max(double(numK),double(max(gIX)))); + +% sort traces by index +im = M; +[nLines,nFrames] = size(im); +[~,I] = sort(gIX); +im = im(I,:); + +%% convert imagesc effect into rgb matrix 'RGB' +if dataFR, % is drawing cell-response fluorescence + cmap = gray(64); + im = AutoScaleImage0to1(im); % scaling to min/max 0/1 + minlim = 0; + maxlim = 1; +else % create blue-white-red map, for regression results + cmap = zeros(64,3); + cmap(:,1) = [linspace(0,1,32), linspace(1,1,32)]; + cmap(:,2) = [linspace(0,1,32), linspace(1,0,32)]; + cmap(:,3) = [linspace(1,1,32), linspace(1,0,32)]; + minlim = -1; %min(min(im)); + maxlim = 1; %max(max(im)); +end +RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap + +%% add vertical color code +if exist('clrmap','var'), + if strcmp(clrmap,'jet'), + temp = flipud(jet(numK)); + else % 'hsv' + temp = hsv(round(numK*1.1)); + end +else % 'hsv' + temp = hsv(round(numK*1.1)); +end +cmap2 = vertcat(temp(1:numK,:),[0,0,0]); % extend colormap to include black +bwidth = max(round(nFrames/30),1); + +idx = gIX(I); +ix_div = [find(diff(idx));length(idx)]; + +bars = ones(nLines,bwidth); +if numK>1, + bars(1:ix_div(1),:) = idx(ix_div(1)); + for i = 2:length(ix_div), + % paint color + bars(ix_div(i-1)+1:ix_div(i),:) = idx(ix_div(i)); + end +end +im_bars = reshape(cmap2(bars,:),[size(bars),3]); +% add a white division margin +div = ones(nLines,round(bwidth/2),3); +% put 3 parts together +im = horzcat(RGB,div,im_bars); + +%% plot figure + +[s1,s2,s3] = size(im); + +hold off; +set(h1,'Position',[pos(1),pos(2)+pos(4)*barratio*3,pos(3),pos(4)*(1-4*barratio)]); +if s1<30, + temp = im; + im = ones(30,s2,s3); + im(1:s1,:,:) = temp; +end +image(im); +set(gca, 'box', 'off') + +hold on; + +% plot cluster division lines +plot([0.5,s2+0.5],[0.5,0.5],'k','Linewidth',0.5); +if numK>1, + for i = 1:length(ix_div),% = numK-1, + y = ix_div(i)+0.5; + plot([0.5,s2+0.5],[y,y],'k','Linewidth',0.5); + end +end + +ylabel(['Number of cells: ' num2str(numcell)]); +set(gca,'YTick',[],'XTick',[]); +set(gcf,'color',[1 1 1]); + +colormap gray; + +% write in number label +x = s2*1.003; +y_last = 0; +for i = 1:length(ix_div), + % avoid label crowding + margin = 0.015*s1; + y0 = ix_div(i)+0.5; + y = max(y_last+margin,y0); + if ilow_high(2)) = low_high(2); +im = mat2gray(temp); +end + +function RGB = ImageToRGB(im,cmap,minlim,maxlim) +L = size(cmap,1); +ix = round(interp1(linspace(minlim,maxlim,L),1:L,im,'linear','extrap')); +RGB = reshape(cmap(ix,:),[size(ix) 3]); % Make RGB image from scaled. +end diff --git a/first version/DrawClustersOnMap_LSh.m b/first version/DrawClustersOnMap_LSh.m new file mode 100644 index 0000000..5bbeb48 --- /dev/null +++ b/first version/DrawClustersOnMap_LSh.m @@ -0,0 +1,137 @@ +function [tot_image, dim_totimage] = DrawClustersOnMap_LSh(CInfo,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,clrmap,full) +%% formatting +[s1,s2] = size(cIX); +if s2>s1, + cIX = cIX'; +end +[s1,s2] = size(gIX); +if s2>s1, + gIX = gIX'; +end + +% down-sample +if ~exist('full','var'), + displaymax = 8000; + if length(cIX) > displaymax, + skip = round(length(cIX)/displaymax); + cIX = cIX(1:skip:end,:); + gIX = gIX(1:skip:end,:); + end +end + +% get numK +if exist('numK','var'), + numK = double(max(numK,max(gIX))); +else + numK = double(max(gIX)); +end + +if strcmp(clrmap,'jet'), + temp = flipud(jet(numK)); +else % 'hsv' + temp = hsv(round(numK*1.1)); +end +cmap = temp(1:numK,:); % extend colormap to include black + +anat_YX = anat_yx/4; +anat_YZ = anat_yz/4; +anat_ZX = anat_zx/4; +anat_ZX = flipud(anat_ZX); +dimv_yx = size(anat_YX); +dimv_yz = size(anat_YZ); +dimv_zx = size(anat_ZX); + +k_zres = 20; +anat_yz2=zeros(dimv_yz(1),dimv_yz(2)*k_zres,3); +anat_zx2=zeros(dimv_zx(1)*k_zres,dimv_zx(2),3); +dim_totimage = [dimv_yx(1)+dimv_zx(1)*k_zres+10,dimv_yx(2)+dimv_yz(2)*k_zres+10,3]; +tot_image=ones(dim_totimage); +% tot_image=zeros(dim_totimage); +% tot_image(:,dimv_yx(2)+(1:10),:)=1; + +% find index manipulation vector to darw circle +circle=makeDisk2(7,15); % make mask of filled circle % (7,15) +mask = zeros(dimv_yx(1),dimv_yx(2)); +mask(1:15,1:15) = circle; +ix = find(mask); +cix = sub2ind([dimv_yx(1),dimv_yx(2)],8,8);% 8 +circle_inds = ix - cix; + +yzplane_inds = -5:5; +zxplane_inds = -5*dimv_zx(1):dimv_zx(1):5*dimv_zx(1); + +weight = 0.3-min(length(cIX)/1000/100,0.1); +for j=1:length(cIX) + if ~isempty(CInfo(cIX(j)).center), + %% Y-X + % cinds = sub2ind([dim_y,dim_x],cell_info(cellsIX(j)).center(1),cell_info(cellsIX(j)).center(2)); + cinds=(CInfo(cIX(j)).center(2)-1)*dimv_yx(1)+CInfo(cIX(j)).center(1); % faster equivalent, lin px idx + labelinds=find((cinds+circle_inds)>0 & (cinds+circle_inds)<=dimv_yx(1)*dimv_yx(2)); % within bounds + ix = gIX(j);%find(U==gIX(j));U = 1:numK; + ixs = cinds+circle_inds(labelinds); + anat_YX(ixs) = cmap(ix,1)*weight + anat_YX(ixs)*(1-weight); % R + ixs = cinds+circle_inds(labelinds)+dimv_yx(1)*dimv_yx(2); + anat_YX(ixs) = cmap(ix,2)*weight + anat_YX(ixs)*(1-weight); % G + ixs = cinds+circle_inds(labelinds)+dimv_yx(1)*dimv_yx(2)*2; + anat_YX(ixs) = cmap(ix,3)*weight + anat_YX(ixs)*(1-weight); % B + + % Y-Z + zweight = weight/2; + % cinds = sub2ind([dim_y,dim_z],cell_info(cellsIX(j)).center(1),cell_info(cellsIX(j)).slice); + cinds=(CInfo(cIX(j)).slice-1)*dimv_yz(1)+CInfo(cIX(j)).center(1); + labelinds=find((cinds+yzplane_inds)>0 & (cinds+yzplane_inds)<=dimv_yz(1)*dimv_yz(2)); + ixs = cinds+yzplane_inds(labelinds); + anat_YZ(ixs) = cmap(ix,1)*zweight + anat_YZ(ixs)*(1-zweight); + ixs = cinds+yzplane_inds(labelinds)+dimv_yz(1)*dimv_yz(2); + anat_YZ(ixs) = cmap(ix,2)*zweight + anat_YZ(ixs)*(1-zweight); + ixs = cinds+yzplane_inds(labelinds)+dimv_yz(1)*dimv_yz(2)*2; + anat_YZ(ixs) = cmap(ix,3)*zweight + anat_YZ(ixs)*(1-zweight); + + % Z-X + zweight = weight/2; + % cinds = sub2ind([dim_y,dim_z],cell_info(cellsIX(j)).center(1),cell_info(cellsIX(j)).slice); + cinds=(CInfo(cIX(j)).center(2)-1)*dimv_zx(1) +(CInfo(cIX(j)).slice); + labelinds=find((cinds+zxplane_inds)>0 & (cinds+zxplane_inds)<=dimv_zx(1)*dimv_zx(2)); + ixs = cinds+zxplane_inds(labelinds); + anat_ZX(ixs) = cmap(ix,1)*zweight + anat_ZX(ixs)*(1-zweight); + ixs = cinds+zxplane_inds(labelinds)+dimv_zx(1)*dimv_zx(2); + anat_ZX(ixs) = cmap(ix,2)*zweight + anat_ZX(ixs)*(1-zweight); + ixs = cinds+zxplane_inds(labelinds)+dimv_zx(1)*dimv_zx(2)*2; + anat_ZX(ixs) = cmap(ix,3)*zweight + anat_ZX(ixs)*(1-zweight); + end +end +% rescale (low-res) z dimension +for k=1:3 + anat_yz2(:,:,1) = imresize(anat_YZ(:,:,1), [dimv_yz(1), dimv_yz(2)*k_zres],'nearest'); + anat_yz2(:,:,2) = imresize(anat_YZ(:,:,2), [dimv_yz(1), dimv_yz(2)*k_zres],'nearest'); + anat_yz2(:,:,3) = imresize(anat_YZ(:,:,3), [dimv_yz(1), dimv_yz(2)*k_zres],'nearest'); + + anat_zx2(:,:,1) = imresize(anat_ZX(:,:,1), [dimv_zx(1)*k_zres, dimv_zx(2)],'nearest'); + anat_zx2(:,:,2) = imresize(anat_ZX(:,:,2), [dimv_zx(1)*k_zres, dimv_zx(2)],'nearest'); + anat_zx2(:,:,3) = imresize(anat_ZX(:,:,3), [dimv_zx(1)*k_zres, dimv_zx(2)],'nearest'); +end + +tot_image(dimv_zx(1)*k_zres+11:end,1:dimv_yx(2),:) = anat_YX; +tot_image(dimv_zx(1)*k_zres+11:end,dimv_yx(2)+11:end,:) = anat_yz2; +tot_image(1:dimv_zx(1)*k_zres,1:dimv_zx(2),:) = flipud(anat_zx2); + +tot_image(tot_image(:)>1) = 1; +tot_image(tot_image(:)<0) = 0; + +image(tot_image); +axis image;axis off + +end + +function out = makeDisk2(radius, dim) +center=floor(dim/2)+1; +out=zeros(dim); +for x=1:dim + for y=1:dim + if norm([x,y]-[center,center])<=radius + out(x,y)=1; + end + end +end + +end \ No newline at end of file diff --git a/FFT_test.m b/first version/FFT_test.m similarity index 100% rename from FFT_test.m rename to first version/FFT_test.m diff --git a/first version/GUI_FishExplorer.m b/first version/GUI_FishExplorer.m new file mode 100644 index 0000000..a417875 --- /dev/null +++ b/first version/GUI_FishExplorer.m @@ -0,0 +1,3446 @@ +%%% Interactive app for exploratory analysis of calcium imaging data +% (with stimulus, behavior, and anatomy) + +% Input calcium data: 1 trace per cell/ROI, ~50,000 cells per fish +% load collection of cells from multiple fish, or load full data of single fish individually +% main outputs: GUI plots, clusters saved into .mat, export variables to MATLAB workspace + +% Tip: to see the structure of this code, use 'Ctrl' + 'm' + '=' to collapse all cells. +% UI controls are organized by tabs and then by rows, instructions and +% comments are where they are constructed ('User Interface:' -> function hfig... ->) +% General internal functions are at the end, some specialized .m functions are outside. + +% Written in Matlab R2014b running on Windows 7. + +% - Xiuye Chen (xiuyechen@gmail.com), Engert Lab, Spring 2015 + + +% To start, run the following in separate script or command window: +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% clear all;close all;clc +% global VAR; +% load('VAR_current.mat','VAR'); +% load('CONSTs_current.mat','CONSTs'); +% hfig = GUI_FishExplorer(CONSTs,VAR); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Full datasets are stored as 'CONST_F?.mat' (? = fish number), to load from the GUI + +%% User Interface: +% function hfig = GUI_FishExplorer(CONSTs,VAR,CONST,i_fish) +function [hfig,fcns] = GUI_FishExplorer(data_dir,CONSTs,VAR,flag_script,var_script) +if exist('flag_script','var') + if ~exist('var_script','var') + var_script={}; + end + runscript(flag_script,var_script); + return +end + +%% Make figure +scrn = get(0,'Screensize'); +hfig = figure('Position',[scrn(3)*0.1 scrn(4)*0.05 scrn(3)*0.85 scrn(4)*0.86],...% [50 100 1700 900] + 'Name','GUI_LSh','DeleteFcn',@closefigure_Callback); +hold off; axis off + +%% Folder setup +% directory for full fish data (.mat) +setappdata(hfig,'data_dir',data_dir); + +% copy of VAR files will be saved into this subfolder: +currentfolder = pwd; +arcmatfolder = [currentfolder '\arc mat']; +if ~exist(arcmatfolder, 'dir') + mkdir(arcmatfolder); +end +setappdata(hfig,'arcmatfolder',arcmatfolder); + +%% Pass external variables into appdata (stored with main figure handle) +setappdata(hfig,'CONSTs',CONSTs); +setappdata(hfig,'VAR',VAR); +nFish = length(CONSTs); + +% parameters / constants +setappdata(hfig,'z_res',19.7); % resoltion in z, um per slice +% fpsec = 1.97; % hard-coded in ext function 'GetStimRegressor.m' +% approx fpsec of 2 used in ext function 'DrawCluster.m' + +% cache +bC = []; % Cache for going b-ack (bad abbr) +fC = []; % Cache for going f-orward +bC.cIX = cell(1,1); +bC.gIX = cell(1,1); +bC.numK = cell(1,1); +fC.cIX = cell(1,1); +fC.gIX = cell(1,1); +fC.numK = cell(1,1); +setappdata(hfig,'bCache',bC); +setappdata(hfig,'fCache',fC); + +% initialization +i_fish = 1; +QuickUpdateFish(hfig,i_fish,'init'); + +%% Initialize internal params into appdata +% fish protocol sets (different sets have different parameters) +M_fish_set = [1, 1, 1, 1, 1, 1, 1, 2, 1]; % M = Matrix +setappdata(hfig,'M_fish_set',M_fish_set); + +% thresholds +thres_merge = 0.9; +thres_split = 0.7; +thres_reg = 0.7; +thres_size = 10; +setappdata(hfig,'thres_merge',thres_merge); +setappdata(hfig,'thres_split',thres_split); % for function 'pushbutton_iter_split' +setappdata(hfig,'thres_reg',thres_reg); % regression threshold, ~correlation coeff +setappdata(hfig,'thres_size',thres_size); % min size for clusters + +% variables +% (not sure all these need to be initialized, probably not complete either) +setappdata(hfig,'clrmap','hsv'); +setappdata(hfig,'opID',0); +setappdata(hfig,'rankID',0); +setappdata(hfig,'rankscore',[]); +setappdata(hfig,'isCentroid',0); +setappdata(hfig,'isWkmeans',1); % in autoclustering, with/without kmeans +setappdata(hfig,'isflipstim',0); % flag for flip updown +setappdata(hfig,'regchoice',{1,1}); % regressor choice; first cell,1=stim,2=motor,3=centroid +setappdata(hfig,'isfullfish',0); % no if QuickUpdateFish, yes if LoadFullFish +setappdata(hfig,'isPlotCorrHist',0); % option for regression +setappdata(hfig,'dataFR',1); % left plot data type, 1 for fluo trace and 0 for regression result +setappdata(hfig,'isPlotReg',1); % plot regressor when selecting it +setappdata(hfig,'hierinplace',1); % hier. partitioning, no reordering + + +Class = getappdata(hfig,'Class'); +classID = numel(Class); +setappdata(hfig,'classID',classID); +setappdata(hfig,'classheader','Test:'); +setappdata(hfig,'newclassname',''); + +Cluster = getappdata(hfig,'Cluster'); +clusID = 1; %numel(Cluster); +setappdata(hfig,'clusID',clusID); +setappdata(hfig,'clusheader','Test:'); +setappdata(hfig,'newclusname',''); + +setappdata(hfig,'cIX',Cluster(clusID).cIX); % cell-IndeX +setappdata(hfig,'gIX',Cluster(clusID).gIX); % group-IndeX + +%% Create UI controls +set(gcf,'DefaultUicontrolUnits','normalized'); +set(gcf,'defaultUicontrolBackgroundColor',[1 1 1]); + +% tab group setup +tgroup = uitabgroup('Parent', hfig, 'Position', [0.05,0.88,0.91,0.12]); +numtabs = 5; +tab = cell(1,numtabs); +M_names = {'General','Operations','Regression','Clustering etc.','Saved Clusters'}; +for i = 1:numtabs, + tab{i} = uitab('Parent', tgroup, 'BackgroundColor', [1,1,1], 'Title', M_names{i}); +end + +% grid setup, to help align display elements +rheight = 0.2; +yrow = 0.7:-0.33:0;%0.97:-0.03:0.88; +bwidth = 0.03; +grid = 0:bwidth+0.001:1; + +% various UI element handles ('global' easier than passing around..) +global hback hfwd hclassname hclassmenu hclusgroupmenu hclusmenu hclusname... + hdatamenu hopID hcentroid hdataFR hloadfish hfishnum hstimreg hmotorreg... + hcentroidreg; + +%% UI ----- tab one ----- (General) +i_tab = 1; + +%% UI row 1: File +i_row = 1; +i = 1;n = 0; + +i=i+n; +n=2; % saves 'VAR' to workspace. 'VAR' contains clustering indices (Class and Clusgroup) for all fish +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Quick save',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_save_Callback); + +i=i+n; +n=2; % saves 'VAR' both to workspace and to 'VAR_current.mat' and to arc folder +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Save .mat',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_savemat_Callback); + +i=i+n; +n=2; % plots selected cells on anatomy z-stack, display and save tiff stack in current directory +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Save Zstack',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_writeZstack_Callback); + +i=i+n; +n=2; % make new figure without the GUI components, can save manually from default menu +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Popup plot',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_popupplot_Callback); + +i=i+n; +n=3; % export main working variables to workspace, can customize! +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Export to workspace',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_exporttoworkspace_Callback); + +%% UI row 2: Load +i_row = 2; +i = 1;n = 0; + +i=i+n; +n=2; % this design is underused now... Quick-load only depends on CONSTs, +% which is a minimum collection of clusters from all fish, so you can load +% the program without full single-fish data. eventually can use this +% platform to do things across fish (like after anatomical alignment). +uicontrol('Parent',tab{i_tab},'Style','text','String','Quick-load fish:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; % loads 'CONSTs_current.mat' from current directory +temp = {}; for j = 1:nFish, temp = [temp,{num2str(j)}];end +hfishnum = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',temp,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@popup_fishmenu_Callback); + +i=i+n; +n=2; % loads full single-fish data from CONST_F?.mat +uicontrol('Parent',tab{i_tab},'Style','text','String','Load full fish:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=2; % these CONST_F?.mat are now saved as uncompressed version 6 .mat, loads faster +temp = {}; for j = 1:nFish, temp = [temp,{num2str(j)}];end +temp = [{'(choose)'},temp]; +hloadfish = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',temp,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@popup_loadfullfishmenu_Callback); + +i=i+n; +n=2; % load different presentations (average, all reps etc), pre-stored in CONST(s) +uicontrol('Parent',tab{i_tab},'Style','text','String','Data type:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=3; +hdatamenu = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',CONSTs{i_fish}.datanames,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@popup_datamenu_Callback); + +i=i+n; +n=4; % only centroids (~mean) of clusters shown on left-side plot, the rest is unchanged +hcentroid = uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Display centroids of clusters',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@checkbox_showcentroids_Callback); + +%% UI ----- tab two ----- (Operations) +i_tab = 2; + +%% UI row 1: Range +i_row = 1; +i = 1;n = 0; + +i=i+n; +n=2; % saves up to 20 steps backwards (datatype/datamenu change does not count) +hback = uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Back',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_back_Callback); + +i=i+n; +n=2; % same, 20 steps forward if applicable +hfwd = uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Forward',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_forward_Callback); + +i=i+n; +n=3; % Choose range of clusters to keep. format: e.g. '1:2,4-6,8:end' +uicontrol('Parent',tab{i_tab},'Style','text','String','Choose cluster range:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_choose_range_Callback); + +i=i+n; +n=2; % Choose range of clusters to exclude. format: e.g. '1:2,4-6,8:end' +uicontrol('Parent',tab{i_tab},'Style','text','String','Exclude:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_exclude_range_Callback); + +i=i+n; +n=1; % Choose range of clusters to fuse/combine into single cluster. format: e.g. '1:2,4-6,8:end' +uicontrol('Parent',tab{i_tab},'Style','text','String','Fuse:',... % (eg 1:2,3-5) + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_fuse_range_Callback); + +i=i+n; +n=2; % can choose range in time, but can't see time axes to decide... resets when choosing datamenu +uicontrol('Parent',tab{i_tab},'Style','text','String','Time-range:',... % (eg 1:2,3-5) + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_t_range_Callback); + +%% UI row 2: Operations +i_row = 2; +i = 1;n = 0; + +i=i+n; +n=2; % operates between the current cell selection and the next (in this order). +uicontrol('Parent',tab{i_tab},'Style','text','String','Set operations:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; % 'setdiff' is current minus next, 'rev setdiff' is next minus current. +n=2; % smartUnion = SmartUnique, cells belonging to 2 clusters goes to the more correlated one +menu = {'(choose)','union','intersect','setdiff','rev setdiff','setxor','smartUnion'}; +hopID = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',menu,'Value',1,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@popup_operations_Callback}); + +i=i+n; +n=2; % rank clusters based on various criteria (to choose) +uicontrol('Parent',tab{i_tab},'Style','text','String','Rank by:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; % 'hier' is the same as default (used after every k-means);'stim-lock' uses std across reps; +n=2; % motor stuff uses the best alignment (by cross-correlation) with the fictive trace; +% L+R is average of L & R; stim-motor is combines 'stim-lock' w 'motor' with arbituary weighting. +menu = {'(choose)','hier.','size','stim-lock','corr',... + 'motor','L motor','R motor','L+R motor','stim-motor'}; +uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',menu,'Value',1,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@popup_ranking_Callback}); + +i=i+n; +n=3; % cluster indices will rank from 1 to number of clusters +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Sqeeze clusters',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_sqeeze_Callback); + +i=i+n; +n=3; % flip the sequenc of clusters, first becomes last +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Flip up-down',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_flipud_Callback); + +i=i+n; +n=3; % switch between 2 colormaps now, jet and a cropped version of hsv (so not all circular) +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Switch colormap',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_clrmap_Callback}); + +%% UI row 3: Anatomy +i_row = 3; +i = 1;n = 0; + +i=i+n; +n=4; % Draw a polygon on anatomy maps to select the cells within those boundaries +uicontrol('Parent',tab{i_tab},'Style','text','String','Draw on anatomy map to crop:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; % Draw on the yx-view (main view) +n=2; % Click to make new vertex, double click to connect to first vertex, +% then optionally drag vertices to reposition, and finally double click again to set +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Draw yx',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_polygon_yx_Callback); + +i=i+n; +n=2; % Draw on yz-view (side projection), same as above +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Draw yz',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_polygon_yz_Callback); + +i=i+n; +n=2; % Draw on zx-view (front projection), same as above +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Draw zx',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_polygon_zx_Callback); + +i=i+n; +n=5; +uicontrol('Parent',tab{i_tab},'Style','text','String','Select all cells within boundaries:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=2; % selects ALL cells contained in dataset that are within the convex shape defined by current cells +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Convex hull',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_withinConvexHull_Callback); + +%% UI ----- tab three ----- (Regression) +i_tab = 3; + +%% UI row 1: regressor +i_row = 1; % Step 1: +i = 1;n = 0; % Choose one type of regressor here, choice highlighted in yellow + +i=i+n; +n=2; % stimulus regressors, go to 'GetStimRegressor.m' to add/update +uicontrol('Parent',tab{i_tab},'Style','text','String','Stim reg.:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=3; % (updated when loading fish) +menu = {'(choose)',''}; +hstimreg = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',menu,'Value',1,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@popup_getstimreg_Callback); + +i=i+n; % motor regressors from fictive, not yet convolved/adjusted for time lag +n=2; % go to 'GetMotorRegressor.m' to add/update +uicontrol('Parent',tab{i_tab},'Style','text','String','Motor reg.:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=2; % (unlike stim regressors, names hardcoded, not importet from regressor...) +menu = {'(choose)','right swims','left swims','forward swims','raw right','raw left','raw average'}; +hmotorreg = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',menu,'Value',1,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@popup_getmotorreg_Callback); + +i=i+n; +n=2; % if checked, plot regressor during selection, together with stim and motor +uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Plot regressor',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],'Value',1,... + 'Callback',@checkbox_popupplotreg_Callback); + +i=i+n+1; +n=4; % Get centroid (~mean) of selected cluster as regressor +uicontrol('Parent',tab{i_tab},'Style','text','String','Regressor from centroid #:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +hcentroidreg = uicontrol('Parent',tab{i_tab},'Style','edit','String',num2str(1),... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_ctrdID_as_reg_Callback); + +i=i+n; +n=2; % too complicated... chomp up centroid, then order as if stimulus were left/right inverted +uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Flip stim',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@checkbox_flipstim_Callback); + +%% UI row 2: regression +i_row = 2; % Step 2: +i = 1;n = 0; % Choose regression, using the regressor chosen above, search in full dataset + +i=i+n; +n=3; +uicontrol('Parent',tab{i_tab},'Style','text','String','Choose regression ->',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=2; % do regression, show all cells with correlation coeff (with regressor) above threshold +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Corr. threshold:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_thres_regression_Callback); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit','String',num2str(thres_reg),... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_regthres_Callback); + +i=i+n; +n=2; % optionally plot histogram of correlation values for all cells in dataset, visualize cut-off +hdataFR = uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Plot corr. hist',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@checkbox_plotcorrhist_Callback); + +i=i+n; +n=3; % probably not so useful. Show top n cells of highest correlation coeff with regressor +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','top corr., number limit:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_topnum_regression_Callback}); + +i=i+n; +n=1; % specify n for above +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_topCorrNumber_Callback); + +i=i+n+1; % more automatic, do a regression with every centroid, then combine (with 'SmartUnique' +n=4; % i.e. overlapping cells are assigned to the cluster with which the correlation is the highest) +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Regression with all centroids',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_AllCentroidRegression_Callback}); + +i=i+n; % this is a remnant button from a failed experiment, idea was to iterate the regression process +n=2; % until the cluster converges, but most of the time it doesn't... +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','iter.reg','Enable','off',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_IterCentroidRegression_Callback}); + +%% UI row 3: (TBD) +i_row = 3; % old idea: display correlation values of multiple regressors (instead of fluo. trace) +i = 1;n = 0; % and then can cluster that and further manipulate... + +i=i+n; +n=4; % did not implement combination of regressors in this version... but display option is coded (dataFR==0) +uicontrol('Parent',tab{i_tab},'Style','text','String','regressor combos??(TBD)',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=4; % also data storage is coded, 'M' could be loaded from M_0_fluo or M_0_reg (storing reg corr values) +hdataFR = uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Display fluo./regression',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@checkbox_dataFR_Callback); + +%% UI ----- tab four ----- (Clustering etc.) +i_tab = 4; + +%% UI row 1: k-means +i_row = 1; +i = 1;n = 0; + +i=i+n; +n=2; % k-means clustering +uicontrol('Parent',tab{i_tab},'Style','text','String','k-means, k =',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_kmeans_Callback); + +i=i+n; +n=4; % anatomy is added to the fluo trace as new dimensions, and (arbituarily) weighted strongly +uicontrol('Parent',tab{i_tab},'Style','text','String','k-means with anatomy, k =',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_kmeans2_Callback); + +i=i+n; % trying to use Silhouette to evaluate cluster quality, find peak to determine optimal k, +n=3; % then display results with that k. But have not set k-means to replicate (speed concern), can be very noisy +uicontrol('Parent',tab{i_tab},'Style','text','String','Find best k in range:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]) + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_kmeans_elbow_Callback); + +%% UI row 2: Auto-clustering +i_row = 2; +i = 1;n = 0; + +i=i+n; % Adjacent clusters (arranged by hier.) will be merged +n=4; % if correlation between centroids is above merging threshold +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Merge thres. (corr. based)',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_merge_Callback}); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit','String',num2str(thres_merge),... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@edit_mergethres_Callback}); + +i=i+n; % further split initial clusters so that the average within-cluster corr coeff is above thres +n=2; % (not so iterative anymore, but could easily restore) +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Iter. split',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_iter_split}); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit','String',num2str(thres_split),... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@edit_splitthres_Callback}); + +i=i+n; +n=3; % minimal size of cluster, otherwise delete at the end +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Cluster size thres.',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_thressize_Callback}); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit','String',num2str(thres_size),... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@edit_sizethres_Callback}); + +i=i+n+1; % longest script here. Splits clusters and prunes them, to yield only very tight clusters. +n=3; % really pretty results, but takes a while when regressing with every centroid. Read code for details. +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Full Auto-Clustering',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_autoclus_Callback}); + +i=i+n; +n=4; % by default it starts with a k-mean of 20 of the current cells. Could skip that if already clustered. +uicontrol('Parent',tab{i_tab},'Style','checkbox','String','(starting with k-mean)',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],'Value',1,... + 'Callback',@checkbox_wkmeans_Callback); + +%% UI row 3: misc plots +i_row = 3; +i = 1;n = 0; + +i=i+n; % k-means are always ranked like in hierachical clustering ~ optimal leaf order +n=2; % here you can rank them again and plot the dengrogram. +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Hier. plot',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_hierplot_Callback); + +i=i+n; % partitioning based on hier. clustering +n=2; % choose between max cluster numbers... +uicontrol('Parent',tab{i_tab},'Style','text','String','Hier.cut, max n:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_hierpartn_Callback); + +i=i+n; % partitioning based on hier. clustering +n=3; % ...or set correlation value threshold +uicontrol('Parent',tab{i_tab},'Style','text','String','Hier.cut, corr thres:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_hierpartthres_Callback); + +i=i+n; % hier. partition in place, i.e. without rearranging order clusters +n=3; +uicontrol('Parent',tab{i_tab},'Style','checkbox','String','Hier.cut in place',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],'Value',1,... + 'Callback',@checkbox_hierinplace_Callback); + +i=i+n; +n=2; % Plots the correlation between all current clusters as a matrix +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Corr. plot',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_corrplot_Callback}); + +i=i+n; +n=2; % random plots not really worth consolidating, only for 1-click convenience +uicontrol('Parent',tab{i_tab},'Style','text','String','temp. plots:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]) + +i=i+n; +n=3; % e.g. this one looks at the history for the fish presented with 4 B/W stimuli +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','plot 4x4(stim)',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_plot4x4_Callback}); + +%% UI ----- tab five ----- (Saved Clusters) +i_tab = 5; + +%% UI row 1: Class +i_row = 1; % There is not really a difference between Class and Cluster +i = 1;n = 0; % I just like to have 2 layers, so I save bigger sets as Class, smaller things as Cluster +% also now Clusters come in cluster-groups ('Clusgroup'), i.e. multiple groups of Clusters... +% the saving is sort of a pain though, several Update___ internal functions dedicated to this. + +i=i+n; +n=1; % +uicontrol('Parent',tab{i_tab},'Style','text','String','Class:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=4; +temp = MakeMenu({Class.name}); +hclassmenu = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',temp,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@popup_classmenu_Callback); + +i=i+n; +n=1; % just edit and press enter to edit current name +uicontrol('Parent',tab{i_tab},'Style','text','String','Edit:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=3; +hclassname = uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@edit_editclassname_Callback}); +set(hclassname,'String',Class(classID).name); + +i=i+n; +n=2; % save the current cells, with current cluster divisions, as 1 class +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Save class',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_saveclass_Callback); + +i=i+n; +n=1; % not really useful, just to help me organize the Class name +uicontrol('Parent',tab{i_tab},'Style','text','String','Header:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit','String','Test',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_class_header_Callback); + +i=i+n; +n=1; % write down new name for new cluster +uicontrol('Parent',tab{i_tab},'Style','text','String','New:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=2; +uicontrol('Parent',tab{i_tab},'Style','edit','String','(blank)',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_newclassname_Callback); + +i=i+n; +n=2; % new class is made when pressing this button +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Make class',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_makeclass_Callback); + +i=i+n; +n=2; % will prompt you to confirm +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Delete Class!',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_delclass_Callback}); + +%% UI row 2: Cluster +i_row = 2; % Virtually the same as Class, except added a ranking option +i = 1;n = 0; + +i=i+n; +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','text','String','Cluster (group):',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +num = length(VAR(i_fish).ClusGroup); +temp = {}; for j = 1:num, temp = [temp,{num2str(j)}];end +hclusgroupmenu = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',temp,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@popup_clusgroupmenu_Callback); + +i=i+n; +n=3; +temp = MakeMenu({Cluster.name}); +hclusmenu = uicontrol('Parent',tab{i_tab},'Style','popupmenu','String',temp,... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@popup_clusmenu_Callback); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','text','String','Edit:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=3; +hclusname = uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@edit_editclusname_Callback}); +set(hclusname,'String',Cluster(clusID).name); + +i=i+n; +n=2; +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Save cluster',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_saveclus_Callback); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','text','String','Header:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit','String','Test',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_clus_header_Callback); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','text','String','New:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=2; +uicontrol('Parent',tab{i_tab},'Style','edit','String','(blank)',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_newclusname_Callback); + +i=i+n; +n=2; +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Make cluster',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@pushbutton_makeclus_Callback); + +i=i+n; +n=2; +uicontrol('Parent',tab{i_tab},'Style','text','String','Set rank:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@edit_setrank_Callback}); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','text','String','Notes:',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=2; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@edit_notes_Callback}); + +i=i+n; +n=2; +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','Delete Cluster!',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_delclus_Callback}); + +%% UI row 3: misc +i_row = 3; % continuation, dealing with the Cluster groups +i = 1;n = 0; % (Cluster-group number: number menu before the Cluster-name menu) + +i=i+n; +n=2; % Combines the chosen clusters into one view (can save as 1 class or cluster then) +uicontrol('Parent',tab{i_tab},'Style','text','String','Union(cluster):',... % (eg 1,3-5) + 'Position',[grid(i) yrow(i_row) bwidth*n rheight]); + +i=i+n; +n=1; +uicontrol('Parent',tab{i_tab},'Style','edit',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',@edit_clusUnion_Callback); + +i=i+n; % just adds a new number to the Clustergroup-number menu, +n=3; % and saves current view as the first cluster in the new Clustergroup +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','new Clustergroup',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_newclusgroup_Callback}); + +i=i+n; +n=3; % delete current +uicontrol('Parent',tab{i_tab},'Style','pushbutton','String','del Clustergroup',... + 'Position',[grid(i) yrow(i_row) bwidth*n rheight],... + 'Callback',{@pushbutton_delclusgroup_Callback}); + +%% Load figure + +UpdateClusID(hfig,clusID); + +%% get local function handles + +fcns = localfunctions; + +end + +%% Callback functions for UI elements: + +%% ----- tab one ----- (General) + +%% row 1: File + +function pushbutton_save_Callback(hObject,~) +hfig = getParentFigure(hObject); +i_fish = getappdata(hfig,'i_fish'); +global VAR; +VAR(i_fish).Class = getappdata(hfig,'Class'); +VAR(i_fish).ClusGroup = CurrentClusGroup(hfig); +disp('saved to workspace ''VAR'''); +end + +function pushbutton_savemat_Callback(hObject,~) +disp('saving...'); +hfig = getParentFigure(hObject); +i_fish = getappdata(hfig,'i_fish'); +arcmatfolder = getappdata(hfig,'arcmatfolder'); +global VAR; +VAR(i_fish).Class = getappdata(hfig,'Class'); +VAR(i_fish).ClusGroup = CurrentClusGroup(hfig); +timestamp = datestr(now,'mmddyy_HHMM'); +matname = [timestamp '.mat']; +save(fullfile(arcmatfolder,matname),'VAR','-v6'); +save('VAR_current.mat','VAR','-v6'); +disp('saved both to workspace and .mat'); +end + +function pushbutton_writeZstack_Callback(hObject,~) +disp('preparing z-stack...'); +isMarkAllCells = 1; + +hfig = getParentFigure(hObject); +isfullfish = getappdata(hfig,'isfullfish'); +if ~isfullfish, + errordlg('Load full fish first!'); + return; +end + +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +CInfo = getappdata(hfig,'CInfo'); +ave_stack = getappdata(hfig,'ave_stack'); + +timestamp = datestr(now,'mmddyy_HHMMSS'); +tiffName = ['stack_' timestamp '.tif']; + +U = unique(gIX); +numK = length(U); + +% left half: stack with cells marked; +% right half: original anatomy, or mark all cells + +ave_stack2=zeros(size(ave_stack,1), size(ave_stack,2)*2, size(ave_stack,3) ,3); +nPlanes=size(ave_stack,3); +dimv_yxz=size(ave_stack); +stacklen=numel(ave_stack); + +circle=makeDisk2(10,21); +[r, v]=find(circle); +r=r-11;v=v-11; +circle_inds = r*dimv_yxz(1)+v; +cmap = hsv(numK); +weight = 0.3; + +for i=1:nPlanes, + ave_stack2(:,:,i,:)=repmat(imNormalize99(ave_stack(:,:,i)),[1 2 1 3]); +end + +for j=1:length(cIX) + cinds=(CInfo(cIX(j)).center(2)-1)*dimv_yxz(1)+CInfo(cIX(j)).center(1); + labelinds=find((cinds+circle_inds)>0 & (cinds+circle_inds)<=dimv_yxz(1)*dimv_yxz(2)); + zinds=dimv_yxz(1)*dimv_yxz(2)*2*(CInfo(cIX(j)).slice-1); + ix = find(U==gIX(j)); + ixs = cinds+circle_inds(labelinds)+zinds; + ave_stack2(ixs)=cmap(ix,1)*weight + ave_stack2(ixs)*(1-weight); + ixs = cinds+circle_inds(labelinds)+zinds+stacklen*2; + ave_stack2(ixs)=cmap(ix,2)*weight + ave_stack2(ixs)*(1-weight); + ixs = cinds+circle_inds(labelinds)+zinds+stacklen*4; + ave_stack2(ixs)=cmap(ix,3)*weight + ave_stack2(ixs)*(1-weight); +end +if isMarkAllCells, + clr = [0,1,0]; + numcell = getappdata(hfig,'numcell'); + for j=1:numcell, + shift = dimv_yxz(1)*dimv_yxz(2); + cinds=(CInfo(j).center(2)-1)*dimv_yxz(1)+CInfo(j).center(1); + labelinds=find((cinds+circle_inds)>0 & (cinds+circle_inds)<=dimv_yxz(1)*dimv_yxz(2)); + zinds=dimv_yxz(1)*dimv_yxz(2)*2*(CInfo(j).slice-1); + ixs = cinds+circle_inds(labelinds)+zinds + shift; + ave_stack2(ixs)=clr(1)*weight + ave_stack2(ixs)*(1-weight); + ixs = cinds+circle_inds(labelinds)+zinds+stacklen*2 + shift; + ave_stack2(ixs)=clr(2)*weight + ave_stack2(ixs)*(1-weight); + ixs = cinds+circle_inds(labelinds)+zinds+stacklen*4 + shift; + ave_stack2(ixs)=clr(3)*weight + ave_stack2(ixs)*(1-weight); + end +end +h = figure; +for i_plane = 1:nPlanes, + im = squeeze(ave_stack2(:,:,i_plane,:)); + image(im); + % save tiff + if (i_plane == 1) + imwrite(im, tiffName, 'compression','none','writemode','overwrite') + else + imwrite(im, tiffName, 'compression','none','writemode','append') + end + pause(0.2) +end +close(h) +end + +function out = makeDisk2(radius, dim) +center=floor(dim/2)+1; +out=zeros(dim); +for x=1:dim + for y=1:dim + if norm([x,y]-[center,center])<=radius + out(x,y)=1; + end + end +end +end + +function pushbutton_popupplot_Callback(hObject,~) +hfig = getParentFigure(hObject); + +% same as function RefreshFigure(hfig), except new axes not global + +% figure('Position',[50,100,853,512]); % right plot size: 2048x1706 +% h1 = axes('Position',[0, 0, 0.5, 1]); % left ~subplot +% h2 = axes('Position',[0.50, 0, 0.5, 1]); % right ~subplot +figure('Position',[50,100,1600,900]); +h1 = axes('Position',[0.05, 0.03, 0.53, 0.94]); % left ~subplot +h2 = axes('Position',[0.61, 0.03, 0.35, 0.94]); % right ~subplot + +CInfo = getappdata(hfig,'CInfo'); +M = getappdata(hfig,'M'); +dataFR = getappdata(hfig,'dataFR'); +fictive = getappdata(hfig,'fictive'); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); + +numK = getappdata(hfig,'numK'); +stim = getappdata(hfig,'stim'); +anat_yx = getappdata(hfig,'anat_yx'); +anat_yz = getappdata(hfig,'anat_yz'); +anat_zx = getappdata(hfig,'anat_zx'); +isCentroid = getappdata(hfig,'isCentroid'); +clrmap = getappdata(hfig,'clrmap'); +rankscore = getappdata(hfig,'rankscore'); +rankID = getappdata(hfig,'rankID'); + +iswrite = (rankID>=2); + +% left subplot +axes(h1); +if isCentroid, + [C,~] = FindCentroid(gIX,M); + DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); +else + DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); +end + +% right subplot +axes(h2); +DrawClustersOnMap_LSh(CInfo,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,clrmap); + +end + +function pushbutton_exporttoworkspace_Callback(hObject,~) +hfig = getParentFigure(hObject); +assignin('base', 'M', getappdata(hfig,'M')); +assignin('base', 'cIX', getappdata(hfig,'cIX')); +assignin('base', 'gIX', getappdata(hfig,'gIX')); +assignin('base', 'numK', getappdata(hfig,'numK')); +assignin('base', 'stim', getappdata(hfig,'stim')); +assignin('base', 'fictive', getappdata(hfig,'fictive')); +assignin('base', 'M_0_fluo', getappdata(hfig,'M_0_fluo')); +end + +%% row 2: Load + +function popup_fishmenu_Callback(hObject,~) +new_i_fish = get(hObject,'Value'); +hfig = getParentFigure(hObject); +QuickUpdateFish(hfig,new_i_fish); +end + +function QuickUpdateFish(hfig,new_i_fish,init) %#ok +CONSTs = getappdata(hfig,'CONSTs'); +setappdata(hfig,'isfullfish',0); + +% load all fields from CONSTs, with names preserved +names = fieldnames(CONSTs{new_i_fish}); % cell of strings +for i = 1:length(names), + setappdata(hfig,names{i},eval(['CONSTs{new_i_fish}.',names{i}])); +end + +% recontruct the 3 big matrices from CIX to original size (numcell) +CIX = CONSTs{new_i_fish}.CIX; +numcell = CONSTs{new_i_fish}.numcell; + +temp = CONSTs{new_i_fish}.CRAZ; +CRAZ = zeros(numcell,size(temp,2)); +CRAZ(CIX,:) = temp; + +temp = CONSTs{new_i_fish}.CRZt; +CRZt = zeros(numcell,size(temp,2)); +CRZt(CIX,:) = temp; + +temp = CONSTs{new_i_fish}.CInfo; +CInfo(numcell).center = ''; +CInfo(numcell).area = ''; +CInfo(numcell).slice = ''; +CInfo(numcell).inds = ''; +CInfo(numcell).y_minmax = ''; +CInfo(numcell).x_minmax = ''; +CInfo(CIX) = temp; + +setappdata(hfig,'CRAZ',CRAZ); % Matrix array, misc collection +setappdata(hfig,'CRZt',CRZt); % Matrix array, misc collection +setappdata(hfig,'CInfo',CInfo); % Matrix array, misc collection + +UpdateFishData(hfig,new_i_fish); +if ~exist('init','var'), + UpdateFishDisplay(hfig,new_i_fish); +end +end + +function UpdateFishData(hfig,new_i_fish) +UpdateDatamenu_Direct(hfig,1); + +% save ClusGroup before updating, if applicable +clusgroupID = getappdata(hfig,'clusgroupID'); +if ~isempty(clusgroupID), + new_clusgroupID = 1; + UpdateClusGroupID(hfig,clusgroupID,new_clusgroupID,'norefresh'); % to save ClusGroup +end +% set new i_fish after UpdateClusGroupID, which saves into old 'i_fish' +setappdata(hfig,'i_fish',new_i_fish); + +% load from VAR +global VAR; +Class = VAR(new_i_fish).Class; % actually collection of classes, each class can have many groups +setappdata(hfig,'Class',Class); + +ClusGroup = VAR(new_i_fish).ClusGroup; % group of 'Cluster', i.e. 2nd level collections +setappdata(hfig,'ClusGroup',ClusGroup); + +clusgroupID = 1; +setappdata(hfig,'clusgroupID',clusgroupID); + +Cluster = ClusGroup{clusgroupID};% collection of clusters, 'Cluster' structurally same as 'Class' +setappdata(hfig,'Cluster',Cluster); +end + +function UpdateFishDisplay(hfig,new_i_fish) +stim = getappdata(hfig,'stim'); +M_fish_set = getappdata(hfig,'M_fish_set'); + +fishset = M_fish_set(new_i_fish); +[~, names] = GetStimRegressor(stim,fishset); +global hstimreg; +set(hstimreg,'String',['(choose)',(names)]); + +% update display +UpdateClassID(hfig,1,'norefresh'); +clusgroupID = 1; +new_clusgroupID = 1; +UpdateClusGroupID(hfig,clusgroupID,new_clusgroupID); % to display new menu +end + +function popup_loadfullfishmenu_Callback(hObject,~) +new_i_fish = get(hObject,'Value')-1; +if new_i_fish>0, + hfig = getParentFigure(hObject); + LoadFullFish(hfig,new_i_fish); +end +global hfishnum; +set(hfishnum,'Value',new_i_fish); +end + +function LoadFullFish(hfig,new_i_fish,CONST) +data_dir=getappdata(hfig,'data_dir'); +if ~exist('CONST','var'), + disp(['loading fish #' num2str(new_i_fish) '...']); + tic + load(fullfile(data_dir,['CONST_F' num2str(new_i_fish) '.mat']),'CONST'); + toc +end +setappdata(hfig,'isfullfish',1); + +% load all fields from CONST, with names preserved +names = fieldnames(CONST); % cell of strings +for i = 1:length(names), + setappdata(hfig,names{i},eval(['CONST.',names{i}])); +end +setappdata(hfig,'numcell',length(CONST.CInfo)); + +UpdateFishData(hfig,new_i_fish); +UpdateFishDisplay(hfig,new_i_fish); +end + +function popup_datamenu_Callback(hObject,~) +ID = get(hObject,'Value'); +hfig = getParentFigure(hObject); +UpdateDatamenu_Direct(hfig,ID); +SetDataFR(hfig,1); % also set M (without updating cIX gIX) +RefreshFigure(hfig); +end + +function UpdateDatamenu_Direct(hfig,ID) +CRAZ = getappdata(hfig,'CRAZ'); +CRZt = getappdata(hfig,'CRZt'); +FcAvr = getappdata(hfig,'FcAvr'); +Fc = getappdata(hfig,'Fc'); +datanames = getappdata(hfig,'datanames'); +tlists = getappdata(hfig,'tlists'); +photostate = getappdata(hfig,'photostate'); + +%% SWITCH LEFT/RIGHT DIRECTION AGAIN +FcAvr = vertcat(FcAvr(2,:),FcAvr(1,:),FcAvr(3:end,:)); +Fc = vertcat(Fc(2,:),Fc(1,:),Fc(3:end,:)); +% now 1=left, 2=right, 3=forward, 4 = raw left, 5 = raw right + +%% +if ID == 1, + M_0 = CRAZ; + fictive = FcAvr; +elseif ID == 2, + M_0 = CRZt; + fictive = Fc; +else + IX = tlists{ID}; + M_0 = CRZt(:,IX); + temp = Fc(:,IX); + % normalize fictive channels, in 2 sets + for k = 1:3,%size(F,1), + m = temp(k,:); + temp(k,:) = (m-min(m))/(max(m)-min(m)); + end + m = temp(4:5,:); + temp(4:5,:) = (m-min(min(m)))/(max(max(m))-min(min(m))); + fictive = temp; +end +% stim from photostate +stim = photostate(tlists{ID}); + +setappdata(hfig,'M_0_fluo',M_0); +setappdata(hfig,'fictive',fictive); +setappdata(hfig,'fictive_0',fictive); +setappdata(hfig,'dataname',datanames(ID)); +setappdata(hfig,'stim',stim); + +global hdatamenu; +if ~isempty(hdatamenu), % before GUI initialization + % 'global' remnant from previous run: hdatamenu is not empty but is invalid object + if isvalid(hdatamenu), + set(hdatamenu,'String',datanames); + set(hdatamenu,'Value',ID); + end +end +end + +function checkbox_showcentroids_Callback(hObject,~) +hfig = getParentFigure(hObject); +setappdata(hfig,'isCentroid',get(hObject,'Value')); +RefreshFigure(hfig); +end + +%% ----- tab two ----- (Operations) + +%% row 1: Range + +function pushbutton_back_Callback(hObject,~) +global hback hfwd; +hfig = getParentFigure(hObject); +bC = getappdata(hfig,'bCache'); +fC = getappdata(hfig,'fCache'); + +if ~isempty(bC.cIX{1}), + % set last into forward cache + fC.cIX = [getappdata(hfig,'cIX'),fC.cIX]; + fC.gIX = [getappdata(hfig,'gIX'),fC.gIX]; + fC.numK = [getappdata(hfig,'numK'),fC.numK]; + % retrieve + cIX = bC.cIX{1}; + gIX = bC.gIX{1}; + numK = bC.numK{1}; + bC.cIX(1) = []; + bC.gIX(1) = []; + bC.numK(1) = []; + % set M + M_0 = GetM_0(hfig); + M = M_0(cIX,:); + % save + setappdata(hfig,'M',M); + setappdata(hfig,'bCache',bC); + setappdata(hfig,'fCache',fC); + setappdata(hfig,'cIX',cIX); + setappdata(hfig,'gIX',gIX); + setappdata(hfig,'numK',numK); + % handle rankID: >=2 means write numbers as text next to colorbar + setappdata(hfig,'rankID',0); + % finish + disp('back (from cache)') + RefreshFigure(hfig); + set(hfwd,'enable','on'); +else % nothing to retrieve + set(hback,'enable','off'); +end +end + +function pushbutton_forward_Callback(hObject,~) +global hback hfwd; +hfig = getParentFigure(hObject); +bC = getappdata(hfig,'bCache'); +fC = getappdata(hfig,'fCache'); + +if ~isempty(fC.cIX{1}), + % set last into (backward) cache + bC.cIX = [getappdata(hfig,'cIX'),bC.cIX]; + bC.gIX = [getappdata(hfig,'gIX'),bC.gIX]; + bC.numK = [getappdata(hfig,'numK'),bC.numK]; + % retrieve + cIX = fC.cIX{1}; + gIX = fC.gIX{1}; + numK = fC.numK{1}; + fC.cIX(1) = []; + fC.gIX(1) = []; + fC.numK(1) = []; + % set M + M_0 = GetM_0(hfig); + M = M_0(cIX,:); + % save + setappdata(hfig,'M',M); + setappdata(hfig,'bCache',bC); + setappdata(hfig,'fCache',fC); + setappdata(hfig,'cIX',cIX); + setappdata(hfig,'gIX',gIX); + setappdata(hfig,'numK',numK); + % handle rankID: >=2 means write numbers as text next to colorbar + setappdata(hfig,'rankID',0); + % finish + disp('forward (from cache)') + RefreshFigure(hfig); + set(hback,'enable','on'); +else + set(hfwd,'enable','off'); +end +end + +function edit_choose_range_Callback(hObject,~) +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +% get/format range +str = get(hObject,'String'); +if ~isempty(str), + str = strrep(str,'end',num2str(max(gIX))); + range = ParseRange(str); + % update indices + tempI = []; + for i = range, + tempI = [tempI;find(gIX==i)]; + end + cIX = cIX(tempI); + gIX = gIX(tempI); + UpdateIndices(hfig,cIX,gIX); + RefreshFigure(hfig); +end +end + +function range = ParseRange(str) +str = strrep(str,':','-'); % e.g. str= '1,3,5:8'; +C = textscan(str,'%d','delimiter',','); +m = C{:}; +range = []; +for i = 1:length(m), + if m(i)>0, + range = [range,m(i)]; %#ok + else % have '-'sign, + range = [range,m(i-1)+1:-m(i)]; %#ok + end +end +end + +function edit_exclude_range_Callback(hObject,~) +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +% get/format range +str = get(hObject,'String'); +if ~isempty(str), + str = strrep(str,'end',num2str(max(gIX))); + range = ParseRange(str); + range = setdiff(unique(gIX),range)'; + % update indices + tempI = []; + for i = range, + tempI = [tempI;find(gIX==i)]; + end + cIX = cIX(tempI); + gIX = gIX(tempI); + UpdateIndices(hfig,cIX,gIX); + RefreshFigure(hfig); +end +end + +function edit_fuse_range_Callback(hObject,~) +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +% get/format range +str = get(hObject,'String'); +if ~isempty(str), + str = strrep(str,'end',num2str(max(gIX))); + range = ParseRange(str); + % update indices + tempI = []; + for i = range, + tempI = [tempI;find(gIX==i)]; + end + gIX(tempI) = gIX(tempI(1)); + [gIX, numK] = SqueezeGroupIX(gIX); + UpdateIndices(hfig,cIX,gIX,numK); + RefreshFigure(hfig); +end +end + +function edit_t_range_Callback(hObject,~) +hfig = getParentFigure(hObject); +M = getappdata(hfig,'M'); + +% get/format range +str = get(hObject,'String'); +if ~isempty(str), + str = strrep(str,'end',num2str(size(M,2))); + range = ParseRange(str); + + UpdateTRange(hfig,range) + RefreshFigure(hfig); +end +end + +function UpdateTRange(hfig,range) +M_0 = GetM_0(hfig); +cIX = getappdata(hfig,'cIX'); +fictive_0 = getappdata(hfig,'fictive_0'); +% photostate_0 = getappdata(hfig,'photostate_0'); + +% set M +M = M_0(cIX,range); +setappdata(hfig,'M',M); +% set fictive +fictive = fictive_0(:,range); +setappdata(hfig,'fictive',fictive); +% % set photostate +% photostate = photostate_0(:,range); +% setappdata(hfig,'photostate',photostate); + +end + +%% row 2: Operations + +function popup_operations_Callback(hObject,~) +opID = get(hObject,'Value') - 1; +hfig = getParentFigure(hObject); +setappdata(hfig,'opID',opID); +% highlight UI +global hopID; +set(hopID,'BackgroundColor',[1,1,0.8]); +end + +function popup_ranking_Callback(hObject,~) +% menu = {'(ranking)','hier.','size','stim-lock','corr','noise'}; +rankID = get(hObject,'Value') - 1; +hfig = getParentFigure(hObject); +setappdata(hfig,'rankID',rankID); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +M = getappdata(hfig,'M'); + +[gIX, numU] = SqueezeGroupIX(gIX); +switch rankID, + case 1, + disp('hier. (default)'); + [gIX, numU] = HierClus(M,gIX); + case 2, + disp('size'); + H = zeros(numU,1); + for i = 1:numU, + H(i) = length(find(gIX==i)); + end + [gIX,rankscore] = SortH(H,gIX,numU,'descend'); + case 3, + disp('stim-lock'); + [gIX,rankscore] = RankByStimLock_Direct(hfig,cIX,gIX,M,numU); + case 4, + disp('corr'); + [~,D] = FindCentroid(gIX,M); + [gIX,rankscore] = SortH(D,gIX,numU); + rankscore = 1-rankscore; + case 5, + disp('motor'); + [gIX,rankscore] = RankByMotorStim_Direct(hfig,gIX,M,numU,1); + case 6, + disp('motor'); + [gIX,rankscore] = RankByMotorStim_Direct(hfig,gIX,M,numU,2); + case 7, + disp('motor'); + [gIX,rankscore] = RankByMotorStim_Direct(hfig,gIX,M,numU,3); + case 8, + disp('motor'); + [gIX,rankscore] = RankByMotorStim_Direct(hfig,gIX,M,numU,4); + case 9, + disp('stim-motor'); + [~,rankscore1] = RankByStimLock_Direct(hfig,cIX,gIX,M,numU); + [~,rankscore2] = RankByMotorStim_Direct(hfig,gIX,M,numU,1); + rankscore = rankscore1 - rankscore2; +end +if rankID>1, + setappdata(hfig,'rankscore',round(rankscore*100)/100); + setappdata(hfig,'clrmap','jet'); +else + setappdata(hfig,'clrmap','hsv'); +end +UpdateIndices(hfig,cIX,gIX,numU); +RefreshFigure(hfig); +disp('ranking complete'); +end + +function [gIX,rankscore] = RankByStimLock_Direct(hfig,cIX,gIX,M,numU) +periods = getappdata(hfig,'periods'); + +if length(periods)==1, + period = periods{1}; + [C,~] = FindCentroid(gIX,M); +else % i_fish>=8, + tlists = getappdata(hfig,'tlists'); + CRZt = getappdata(hfig,'CRZt'); + IX = tlists{6}; % ptomr_circ + M_ = CRZt(cIX,IX); + [C,~] = FindCentroid(gIX,M_); + periods = getappdata(hfig,'periods'); + period = periods{1}+periods{2}; +end +% C2 = zscore(C,0,2); +% C_3D = reshape(C2,size(C2,1),period,[]); +C_3D_0 = reshape(C,size(C,1),period,[]); +C_3D = zscore(C_3D_0,0,2); + +H = nanmean(nanstd(C_3D,0,3),2); +[gIX,rankscore] = SortH(H,gIX,numU); +end + +function [gIX,rankscore] = RankByMotorStim_Direct(hfig,gIX,M,numU,option) +C = FindCentroid(gIX,M); +fictive = getappdata(hfig,'fictive'); + +reg = zeros(3,length(fictive)); +reg(1,:) = fictive(5,:); % Left +reg(2,:) = fictive(4,:); % Right +reg(3,:) = mean(fictive(4:5,:)); +H = zeros(numU,1); +shift = zeros(numU,1); +a = zeros(1,3); +I = zeros(1,3); +for i = 1:numU, + switch option, + case 1, + for j = 1:3, + [a(j),I(j)] = max(abs(xcorr(C(i,:),reg(j,:),'coeff'))); + end + [H(i),jmax] = max(a); + shift(i) = I(jmax) - length(fictive); + case 2, + [H(i),I] = max(xcorr(C(i,:),reg(1,:),'coeff')); + shift(i) = I - length(fictive); + case 3, + [H(i),I] = max(xcorr(C(i,:),reg(2,:),'coeff')); + shift(i) = I - length(fictive); + case 4, + [H(i),I] = max(xcorr(C(i,:),reg(3,:),'coeff')); + shift(i) = I - length(fictive); + end +end +[gIX,rankscore] = SortH(H,gIX,numU,'descend'); +assignin('base', 'shift', shift); +end + +function [gIX,B] = SortH(H,gIX,numU,descend) % new gIX is sorted based on H, size(H)=[numU,1]; +if exist('descend','var'), + [B,I] = sort(H,'descend'); +else + [B,I] = sort(H); +end +gIX_last = gIX; +for i = 1:numU, + gIX(gIX_last==I(i)) = i; +end +end + +function pushbutton_sqeeze_Callback(hObject,~) +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +[gIX, numU] = SqueezeGroupIX(gIX); +UpdateIndices(hfig,cIX,gIX,numU); +RefreshFigure(hfig); +end + +function pushbutton_flipud_Callback(hObject,~) +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); + +U = unique(gIX); +U_ = flipud(U); + +gIX_ = gIX; +for i = 1:length(U), + gIX_(gIX==U(i)) = U_(i); +end + +UpdateIndices(hfig,cIX,gIX_); +RefreshFigure(hfig); +end + +function pushbutton_clrmap_Callback(hObject,~) +hfig = getParentFigure(hObject); +clrmap = getappdata(hfig,'clrmap'); +if strcmp(clrmap,'jet'), + setappdata(hfig,'clrmap','hsv'); +else + setappdata(hfig,'clrmap','jet'); +end +RefreshFigure(hfig); +end + +%% row 3: Anatomy + +function pushbutton_polygon_yx_Callback(hObject,~) +k_zres = 20; + +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +% M = getappdata(hfig,'M'); +CInfo = getappdata(hfig,'CInfo'); +numcell = getappdata(hfig,'numcell'); +anat_yx = getappdata(hfig,'anat_yx'); +anat_zx = getappdata(hfig,'anat_zx'); +dimv_yx = size(anat_yx); +dimv_zx = size(anat_zx); +isfullfish = getappdata(hfig,'isfullfish'); + +h_poly_yx = impoly; +wait(h_poly_yx); % double click to finalize position! +% update finalized polygon in bright color +setColor(h_poly_yx,[0 1 1]); + +if isfullfish, + IJs = reshape([CInfo.center],2,[])'; +else % Matlab can't handle CInfo with all the empty entries -> manual padding + IJs_ = reshape([CInfo(cIX).center],2,[])'; + IJs = ones(numcell,2); + IJs(cIX,:) = IJs_; +end +A = sub2ind(dimv_yx(1:2),IJs(:,1),IJs(:,2)); +MaskArray = createMask(h_poly_yx); +MaskArray(1:dimv_zx*k_zres+10,:) = []; +B = find(MaskArray); % find indices of pixels within ROI + +temp = find(ismember(A,B)); +[IX,ia,~] = intersect(cIX,temp); + +cIX = IX; +gIX = gIX(ia); + +UpdateIndices(hfig,cIX,gIX); +RefreshFigure(hfig); +end + +function pushbutton_polygon_yz_Callback(hObject,~) +k_zres = 20; + +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +% M = getappdata(hfig,'M'); +CInfo = getappdata(hfig,'CInfo'); +numcell = getappdata(hfig,'numcell'); +anat_yx = getappdata(hfig,'anat_yx'); +anat_yz = getappdata(hfig,'anat_yz'); +anat_zx = getappdata(hfig,'anat_zx'); +dimv_yx = size(anat_yx); +dimv_yz = size(anat_yz); +dimv_zx = size(anat_zx); +isfullfish = getappdata(hfig,'isfullfish'); + +h_poly_z = impoly; +wait(h_poly_z); % double click to finalize position! +% update finalized polygon in bright color +setColor(h_poly_z,[0 1 1]); + +if isfullfish, + IJs = reshape([CInfo.center],2,[])'; + Zs = [CInfo.slice]'; +else % Matlab can't handle CInfo with all the empty entries -> manual padding + IJs_ = reshape([CInfo(cIX).center],2,[])'; + IJs = ones(numcell,2); + IJs(cIX,:) = IJs_; + Zs_ = [CInfo(cIX).slice]'; + Zs = ones(numcell,1); + Zs(cIX) = Zs_; +end +A = sub2ind(dimv_yz(1:2),IJs(:,1),Zs); + +MaskArray = createMask(h_poly_z); +MaskArray(1:dimv_zx*k_zres+10,:) = []; +MaskArray(:,1:dimv_yx(2)+10) = []; +zMaskArray = imresize(MaskArray,[dimv_yz(1),dimv_yz(2)]); +B = find(zMaskArray); % find indices of pixels within ROI + +temp = find(ismember(A,B)); +[IX,ia,~] = intersect(cIX,temp); + +cIX = IX; +gIX = gIX(ia); + +UpdateIndices(hfig,cIX,gIX); +RefreshFigure(hfig); +end + +function pushbutton_polygon_zx_Callback(hObject,~) +k_zres = 20; + +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +% M = getappdata(hfig,'M'); +CInfo = getappdata(hfig,'CInfo'); +numcell = getappdata(hfig,'numcell'); +% anat_yx = getappdata(hfig,'anat_yx'); +anat_zx = getappdata(hfig,'anat_zx'); +% dimv_yx = size(anat_yx); +dimv_zx = size(anat_zx); +isfullfish = getappdata(hfig,'isfullfish'); + +h_poly_z = impoly; +wait(h_poly_z); % double click to finalize position! +% update finalized polygon in bright color +setColor(h_poly_z,[0 1 1]); + +if isfullfish, + IJs = reshape([CInfo.center],2,[])'; + Zs = [CInfo.slice]'; +else % Matlab can't handle CInfo with all the empty entries -> manual padding + IJs_ = reshape([CInfo(cIX).center],2,[])'; + IJs = ones(numcell,2); + IJs(cIX,:) = IJs_; + Zs_ = [CInfo(cIX).slice]'; + Zs = ones(numcell,1); + Zs(cIX) = Zs_; +end +A = sub2ind(dimv_zx(1:2),Zs,IJs(:,2)); + +MaskArray = createMask(h_poly_z); + +MaskArray(dimv_zx(1)*k_zres+1:end,:) = []; +MaskArray(:,dimv_zx(2)+1:end) = []; +zMaskArray = imresize(MaskArray,[dimv_zx(1),dimv_zx(2)]); +zMaskArray = flipud(zMaskArray); % careful! +B = find(zMaskArray); % find indices of pixels within ROI + +temp = find(ismember(A,B)); +[IX,ia,~] = intersect(cIX,temp); + +cIX = IX; +gIX = gIX(ia); + +UpdateIndices(hfig,cIX,gIX); +RefreshFigure(hfig); +end + +function pushbutton_withinConvexHull_Callback(hObject,~) % only works with full load +disp('find points within boundary...'); +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +CInfo = getappdata(hfig,'CInfo'); +IJs = reshape([CInfo.center],2,[])'; +Zs = [CInfo.slice]'; +% point sets +xyz_all = horzcat(IJs, Zs); +xyz_cIX = horzcat(IJs(cIX,:), Zs(cIX)); + +tri = delaunayn(xyz_cIX); % Generate delaunay triangulization +t = tsearchn(xyz_cIX, tri, xyz_all); % Determine which triangle point is within +I_inside = ~isnan(t); + +cIX_1 = cIX; +gIX_1 = gIX; +cIX_2 = setdiff(find(I_inside),cIX_1); +gIX_2 = ones(length(cIX_2),1)*(max(gIX_1)+1); +cIX = [cIX_1;cIX_2]; +gIX = [gIX_1;gIX_2]; + +UpdateIndices(hfig,cIX,gIX,length(unique(gIX))); +RefreshFigure(hfig); +disp('search complete'); +end + +%% ----- tab three ----- (Regression) + +%% row 1: regressor + +function popup_getstimreg_Callback(hObject,~) +i_reg = get(hObject,'Value')-1; +hfig = getParentFigure(hObject); +if i_reg==0, + return; +end +setappdata(hfig,'regchoice',{1,i_reg}); +% highlight the choice (yellow) +global hstimreg hmotorreg hcentroidreg; +set(hstimreg,'BackgroundColor',[1,1,0.8]); % yellow +set(hmotorreg,'BackgroundColor',[1,1,1]); +set(hcentroidreg,'BackgroundColor',[1,1,1]); + +isPlotReg = getappdata(hfig,'isPlotReg'); +if isPlotReg, + PlotRegWithStimMotor(hfig); +end +end + +function PlotRegWithStimMotor(hfig) +stim = getappdata(hfig,'stim'); +fictive = getappdata(hfig,'fictive'); +regressor = GetRegressor(hfig); +% plot regressor +figure; +halfbarheight = 1; +stimbar = GetStimBar(halfbarheight,stim); +reg = imNormalize99(regressor); +regbar = repmat(reg,[1,1,3]); + +subplot(3,1,1);image(stimbar); axis off; title('stimulus'); +subplot(3,1,2);image(regbar); axis off; title('regressor'); +subplot(3,1,3);imagesc(fictive);colormap hot; axis off; title('motor'); +end + +function out=imNormalize99(im) +im=double(im); +temp=sort(im(:),'descend'); +th1=temp(round(length(im(:))/100)); +th2=min(im(:)); + +out=(im-th2)/(th1-th2); +out(out>1)=1; +end + +function popup_getmotorreg_Callback(hObject,~) +hfig = getParentFigure(hObject); +i_reg = get(hObject,'Value')-1; +% {'left swims','right swims','forward swims','raw left','raw right','raw average'}; +if i_reg>0, + setappdata(hfig,'regchoice',{2,i_reg}); +end +% highlight the choice +global hstimreg hmotorreg hcentroidreg; +set(hstimreg,'BackgroundColor',[1,1,1]); +set(hmotorreg,'BackgroundColor',[1,1,0.8]); % yellow +set(hcentroidreg,'BackgroundColor',[1,1,1]); + +isPlotReg = getappdata(hfig,'isPlotReg'); +if isPlotReg, + PlotRegWithStimMotor(hfig); +end +end + +function checkbox_popupplotreg_Callback(hObject,~) +hfig = getParentFigure(hObject); +setappdata(hfig,'isPlotReg',get(hObject,'Value')); +end + +function edit_ctrdID_as_reg_Callback(hObject,~) +str = get(hObject,'String'); +if ~isempty(str), + temp = textscan(str,'%d'); + ctrdID = temp{:}; +end +hfig = getParentFigure(hObject); +setappdata(hfig,'regchoice',{3,ctrdID}); +% highlight the choice +global hstimreg hmotorreg hcentroidreg; +set(hstimreg,'BackgroundColor',[1,1,1]); +set(hmotorreg,'BackgroundColor',[1,1,1]); +set(hcentroidreg,'BackgroundColor',[1,1,0.8]); % yellow +end + +function checkbox_flipstim_Callback(hObject,~) +hfig = getParentFigure(hObject); +setappdata(hfig,'isflipstim',get(hObject,'Value')); +end + +%% row 2: regression + +function regressor = GetRegressor(hObject) +hfig = getParentFigure(hObject); +regchoice = getappdata(hfig,'regchoice'); +stim = getappdata(hfig,'stim'); + +if regchoice{1}==1, % stim Regressor + i_fish = getappdata(hfig,'i_fish'); + M_fish_set = getappdata(hfig,'M_fish_set'); + fishset = M_fish_set(i_fish); + [regressors, names] = GetStimRegressor(stim,fishset); + regressor = regressors(regchoice{2}).im; + +elseif regchoice{1}==2, % motor Regressor + fictive = getappdata(hfig,'fictive'); + regressors = GetMotorRegressor(fictive); + regressor = regressors(regchoice{2}).im; + +else % regchoice{1}==3, from Centroid + ctrdID = regchoice{2}; + isflipstim = getappdata(hfig,'isflipstim'); + M = getappdata(hfig,'M'); + gIX = getappdata(hfig,'gIX'); + + i = find(unique(gIX)==ctrdID); + if isempty(i), + disp('input is empty!');beep; + regressor = []; + return; + end + [C,~] = FindCentroid(gIX,M); + regressor = C(i,:); + + if isflipstim, % swap left/right stim + % 0 = all black; 1 = black/white; 2 = white/black; 3 = all white; 4 = all gray; + % 10 = forward grating (very slow, more for calibration) + % 11 = rightward grating + % 12 = leftward grating + % swap: 1~2,11~12 + + % if i_fish == 6 || i_fish == 7, + % halflength = length(stim)/2; + % stim = stim(1:halflength); % length is always even + % end + flipstim = stim; + flipstim(stim==1)=2; + flipstim(stim==2)=1; + flipstim(stim==11)=12; + flipstim(stim==12)=11; + + phototrans = GetPhotoTrans(stim); + fliptrans = GetPhotoTrans(flipstim); + + % transition e.g.: + % photostate: 2 3 1 3 3 0 0 1 1 0 3 2 0 2 2 1 + % phototrans: 6 11 13 7 15 12 0 1 5 4 3 14 8 2 10 9 + + IX = []; + targetstates = unique(fliptrans,'stable'); + for i = 1:length(targetstates), + IX = [IX, find(phototrans==targetstates(i))]; + end + % if i_fish == 6 || i_fish == 7, + % reg = regressor(1:halflength); + % reg = reg(IX); + % regressor = repmat(reg,1,2); + % else + regressor = regressor(IX); + % end + end +end + +end + +function checkbox_plotcorrhist_Callback(hObject,~) +hfig = getParentFigure(hObject); +setappdata(hfig,'isPlotCorrHist',get(hObject,'Value')); +end + +function pushbutton_topnum_regression_Callback(hObject,~) +disp('regression...'); +hfig = getParentFigure(hObject); +M_0 = getappdata(hfig,'M_0_fluo'); +numTopCorr = getappdata(hfig,'numTopCorr'); +regressor = GetRegressor(hfig); + +%% for each cell, find correlation coeff +M_corr = corr(regressor',M_0');%cell_resp_ave'); +% M_corr_ctrl = corr(regressor.ctrl',M_0'); + +[~,I] = sort(M_corr,'descend'); +cIX = I(1:numTopCorr)'; +gIX = ceil((1:length(cIX))'/length(cIX)*min(20,length(cIX))); + +score = M_corr(cIX); +rankscore = zeros(1,20); +for i = 1:20, + ix = find(gIX==i); + rankscore(i) = mean(score(ix)); +end +setappdata(hfig,'rankscore',round(rankscore*100)/100); +setappdata(hfig,'rankID',99); + +UpdateIndices(hfig,cIX,gIX,20); +RefreshFigure(hfig); +end + +function edit_topCorrNumber_Callback(hObject,~) +str = get(hObject,'String'); +if ~isempty(str), + temp = textscan(str,'%d'); + numTopCorr = temp{:}; + hfig = getParentFigure(hObject); + setappdata(hfig,'numTopCorr',numTopCorr); +end +end + +function pushbutton_thres_regression_Callback(hObject,~) +disp('regression...'); +hfig = getParentFigure(hObject); +M_0 = getappdata(hfig,'M_0_fluo'); +thres_reg = getappdata(hfig,'thres_reg'); +regressor = GetRegressor(hfig); +if isempty(regressor), + return; +end + +isPlotCorrHist = getappdata(hfig,'isPlotCorrHist'); + +[cIX,gIX] = Regression_direct(M_0,thres_reg,regressor,isPlotCorrHist); + +if ~isempty(gIX), + gIX = ceil((1:length(cIX))'/length(cIX)*min(20,length(cIX))); + UpdateIndices(hfig,cIX,gIX,40); + RefreshFigure(hfig); +end +end + +function [cIX_,gIX_,R] = Regression_direct(M_0,thres_reg,regressor,isPlotCorrHist) % gIX_ is just ones +R = corr(regressor',M_0'); +if thres_reg>0, + cIX_ = (find(R>thres_reg))'; +elseif thres_reg<0, + cIX_ = (find(R + cIX_ = Regression_direct(M_0,thres_reg,regressor); + clusters(i).cIX = cIX_; +end +disp('union...'); +CIX = vertcat(clusters(1:end).cIX); +GIX = []; +for i = 1:numel(clusters), + GIX = [GIX; ones(length(clusters(i).cIX),1)*i]; %#ok +end +[cIX,gIX,numK] = SmartUnique(CIX,GIX,M_0(CIX,:)); +end + +function pushbutton_IterCentroidRegression_Callback(hObject,~) +hfig = getParentFigure(hObject); +M_0 = getappdata(hfig,'M_0_fluo'); +thres_reg = getappdata(hfig,'thres_reg'); +M = getappdata(hfig,'M'); +gIX = getappdata(hfig,'gIX'); +regchoice = getappdata(hfig,'regchoice'); +if regchoice{1}==3, % from Centroid + ctrdID = regchoice{2}; +end + +i = find(unique(gIX)==ctrdID); +if isempty(i), + disp('input is empty!');beep; + return; +end +[C,~] = FindCentroid(gIX,M); +regressor = C(i,:); +[cIX_,gIX_] = Regression_direct(M_0,thres_reg,regressor); +numC = length(cIX_); +regressor_last = regressor; + +if ~isempty(cIX_), + for itr = 1:20, + [cIX_,gIX_] = Regression_direct(M_0,thres_reg,regressor_last); + numC = length(cIX_); + disp(num2str(numC)); + M = M_0(cIX_,:); + regressor = FindCentroid(gIX_,M); + if itr>1 && abs(numC-numC_last)<5,%corr(regressor_last',regressor')>thres_iter, + disp(['converged at itr = ' num2str(itr)]); + break; + end + regressor_last = regressor; + numC_last = numC; + end + gIX_ = ceil((1:length(cIX_))'/length(cIX_)*min(20,length(cIX_))); + UpdateIndices(hfig,cIX_,gIX_,40); + RefreshFigure(hfig); +end +end + +%% row 3: (TBD) + +function checkbox_dataFR_Callback(hObject,~) +hfig = getParentFigure(hObject); +SetDataFR(hfig,get(hObject,'Value')); +RefreshFigure(hfig); +end + +%% ----- tab four ----- (Clustering etc.) + +%% row 1: k-means + +function edit_kmeans_Callback(hObject,~) +hfig = getParentFigure(hObject); +M = getappdata(hfig,'M'); + +str = get(hObject,'String'); +if ~isempty(str), + temp = textscan(str,'%d'); + numK = temp{:}; + disp(['k-means k=' num2str(numK) '...']); + tic + rng('default');% default = 0, but can try different seeds if doesn't converge + if numel(M)*numK < 10^7 && numK~=1, + disp('Replicates = 5'); + [gIX,C,sumd,D] = kmeans(M,numK,'distance','correlation','Replicates',5); + elseif numel(M)*numK < 10^8 && numK~=1, + disp('Replicates = 3'); + [gIX,C,sumd,D] = kmeans(M,numK,'distance','correlation','Replicates',3); + else + [gIX,C,sumd,D] = kmeans(M,numK,'distance','correlation');%,'Replicates',3); + end + toc + beep + + if numK>1, + gIX = HierClusDirect(C,gIX,numK); + end + + cIX = getappdata(hfig,'cIX'); + UpdateIndices(hfig,cIX,gIX,numK); + RefreshFigure(hfig); +end +end + +function edit_kmeans2_Callback(hObject,~) % based on anatomical distance +hfig = getParentFigure(hObject); +z_res = getappdata(hfig,'z_res'); +cIX = getappdata(hfig,'cIX'); +CInfo = getappdata(hfig,'CInfo'); +M = getappdata(hfig,'M'); + +XY = reshape([CInfo(cIX).center],2,[])'; +Z = [CInfo(cIX).slice]'.*z_res; +XYZ = horzcat(XY,Z); +Combo = horzcat(XYZ/50,M);% how to weight???????????????? + +str = get(hObject,'String'); +if ~isempty(str), + temp = textscan(str,'%d'); + numK = temp{:}; + disp(['anat. weighted k-means k=' num2str(numK) '...']); + tic + rng('default');% default = 0, but can try different seeds if doesn't converge + [gIX,C,sumd,D] = kmeans(Combo,numK,'distance','correlation');%,'Replicates',3); + toc + beep + + if numK>1, + gIX = HierClusDirect(C,gIX,numK); + end + + cIX = getappdata(hfig,'cIX'); + UpdateIndices(hfig,cIX,gIX,numK); + RefreshFigure(hfig); +end +end + +function edit_kmeans_elbow_Callback(hObject,~) +hfig = getParentFigure(hObject); +gIX = getappdata(hfig,'gIX'); +M = getappdata(hfig,'M'); + +str = get(hObject,'String'); +if ~isempty(str), + str = strrep(str,'end',num2str(max(gIX))); + range = ParseRange(str); + + Silh_mean = zeros(1,length(range)); + for i = 1:length(range), + disp(['k-means k=' num2str(range(i))]); + gIX = kmeans(M,range(i),'distance','correlation'); + silh = silhouette(M,gIX,'correlation'); + Silh_mean(i) = mean(silh); + end + [~,ix] = max(Silh_mean); + numK = range(ix); + + figure;hold on + plot(range,Silh_mean,'o-','color',[0.5 0.5 0.5]); + + if length(range)~=1, + plot([numK,numK],[min(Silh_mean),max(Silh_mean)],'r--'); + xlim([range(1),range(end)]); + ylim([min(Silh_mean),max(Silh_mean)]) + end + xlabel('k = # of clusters') + ylabel('silhouette mean') + + %% perform regular k-means + disp(['k-means k=' num2str(numK)]); + rng('default');% default = 0, but can try different seeds if doesn't converge + if numel(M)*numK < 10^7 && numK~=1, + disp('Replicates = 5'); + [gIX,C,sumd,D] = kmeans(M,numK,'distance','correlation','Replicates',5); + elseif numel(M)*numK < 10^8 && numK~=1, + disp('Replicates = 3'); + [gIX,C,sumd,D] = kmeans(M,numK,'distance','correlation','Replicates',3); + else + [gIX,C,sumd,D] = kmeans(M,numK,'distance','correlation');%,'Replicates',3); + end + if numK>1, + gIX = HierClusDirect(C,gIX,numK); + end + cIX = getappdata(hfig,'cIX'); + UpdateIndices(hfig,cIX,gIX,numK); + RefreshFigure(hfig); + beep; +end +end + +%% row 2: Auto-clustering + +function pushbutton_merge_Callback(hObject,~) +% disp('merging...'); +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +M = getappdata(hfig,'M'); +U = unique(gIX); +numU = length(U); +[C,D] = FindCentroid(gIX,M); + +thres_merge = getappdata(hfig,'thres_merge'); + +i = 1; +while i thres_merge, + IX = find(gIX == U(i+1)); + gIX(IX)=U(i); + U = unique(gIX); + numU = length(U); + + IX = find(gIX == U(i)); + M_s = M(IX,:); + [~,C1,~,D1] = kmeans(M_s,1,'distance','correlation'); + C(i,:) = C1; + D(i) = mean(D1); + C(i+1,:) = []; + D(i+1) = []; + else + i = i+1; + end +end + +if numU>1, + [gIX, numU] = HierClus(M,gIX); +end + +UpdateIndices(hfig,cIX,gIX,numU); +RefreshFigure(hfig); +disp('merging complete'); +end + +function edit_mergethres_Callback(hObject,~) +str = get(hObject,'String'); +temp = textscan(str,'%f',1); +hfig = getParentFigure(hObject); +setappdata(hfig,'thres_merge',temp{:}); +end + +function pushbutton_iter_split(hObject,~) +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +numU = getappdata(hfig,'numK'); +thres_split = getappdata(hfig,'thres_split'); +M_0 = getappdata(hfig,'M_0_fluo'); +classheader = getappdata(hfig,'classheader'); + +disp('iter. split all, beep when done...'); +thres_size = 10; +thres_H = thres_split; +% thres_H = [0.2;0.15;0.1;0.05]; % could have more rounds... + +% initialization +I_rest = []; +% loop +tic +for round = 1:length(thres_H), + disp(['round ' num2str(round) ', numU = ' num2str(numU)]); + dthres = 1-thres_H(round); + gIX_last = gIX; + I_clean_last = cIX; + cIX = []; + gIX = []; + for i = 1:numU, + disp(['i = ' num2str(i)]); + ix = find(gIX_last == i); +% IX = ix; + IX = I_clean_last(ix); + M_s = M_0(IX,:); + [I_rest,cIX,gIX,numU] = CleanClus(M_s,IX,I_rest,cIX,gIX,numU,dthres,thres_size); +% cIX = I_clean_last(I_clean); + end + + [gIX, numU] = SqueezeGroupIX(gIX); + SaveClass(hfig,cIX,gIX,classheader,['clean_round' num2str(round)]); + + SaveClass(hfig,I_rest,ones(length(I_rest),1),classheader,... + ['rest_round' num2str(round)]); +end +toc +beep +end + +function [I_rest,I_clean,gIX_clean,numU] = CleanClus(M_s,IX,I_rest,I_clean,gIX_clean,numU,dthres,thres_size) +I_clean_s = []; + +% find numK_s for kmeans +kmax = min(round(size(M_s,1)/thres_size),30); +% try numK_s = 1 +numK_s = 1; +rng('default'); +[gIX_s,~,~,D] = kmeans(M_s,numK_s,'distance','correlation'); +Dist = min(D,[],2); +if mean(Dist)>dthres, + % try numK_s = kmax + numK_s = kmax; + rng('default'); + [gIX_s,~,~,D] = kmeans(M_s,numK_s,'distance','correlation'); + Dist = min(D,[],2); + if mean(Dist)thres_silh, + % keep the k-means k=2 subsplit + disp('split'); + gIX(IX) = gIX_ + numU; % reassign (much larger) gIX + end +end +[gIX, ~] = SqueezeGroupIX(gIX); + +SaveClass(hfig,cIX,gIX,classheader,'clean_round3'); + +%% rank by stim-lock ?? bug? +% disp('stim-lock'); +% M = M_0(cIX,:); +% [gIX,rankscore] = RankByStimLock_Direct(hfig,cIX,gIX,M,numU); +% disp('ranking complete'); +% % and threshold +% IX = find(rankscore thres_merge, + IX = find(gIX == U(i+1)); + gIX(IX)=U(i); + U = unique(gIX); + numU = length(U); + + IX = find(gIX == U(i)); + M_s = M(IX,:); + [~,C1,~,D1] = kmeans(M_s,1,'distance','correlation'); + C(i,:) = C1; + D(i) = mean(D1); + C(i+1,:) = []; + D(i+1) = []; + else + i = i+1; + end +end +[gIX, numU] = HierClus(M,gIX); +disp('merging complete'); +end + +function checkbox_wkmeans_Callback(hObject,~) +hfig = getParentFigure(hObject); +setappdata(hfig,'isWkmeans',get(hObject,'Value')); +end + +%% row 3: misc plots + +function pushbutton_hierplot_Callback(hObject,~) +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +M = getappdata(hfig,'M'); + +[gIX2, numU] = HierClus(M,gIX,'isplotfig'); +if ~isequal(gIX,gIX2), + UpdateIndices(hfig,cIX,gIX2,numU); + RefreshFigure(hfig); +end +end + +function edit_hierpartn_Callback(hObject,~) +hfig = getParentFigure(hObject); +M = getappdata(hfig,'M'); +gIX = getappdata(hfig,'gIX'); +hierinplace = getappdata(hfig,'hierinplace'); + +[gIX, numU] = SqueezeGroupIX(gIX); +[C,~] = FindCentroid(gIX,M); + +str = get(hObject,'String'); +if ~isempty(str), + temp = textscan(str,'%d',1); + numCuts = temp{:}; + +% tree = linkage(C,'average','correlation'); +% IX_tree = cluster(tree,'maxclust',numCuts); + IX_tree = clusterdata(C,'criterion','distance','distance','correlation','maxclust',numCuts) + + if hierinplace, + % sort to keep clusters in place + U = unique(IX_tree,'stable'); + temp = zeros(size(IX_tree)); + for i = 1:length(U), + temp(IX_tree==U(i)) = i; + end + IX_tree = temp; + end + + % update gIX + temp = zeros(size(gIX)); + for i = 1:numU, + temp(gIX==i) = IX_tree(i); + end + gIX = temp; + + cIX = getappdata(hfig,'cIX'); + UpdateIndices(hfig,cIX,gIX,length(unique(IX_tree))); + RefreshFigure(hfig); +end +end + +function edit_hierpartthres_Callback(hObject,~) +hfig = getParentFigure(hObject); +M = getappdata(hfig,'M'); +gIX = getappdata(hfig,'gIX'); +hierinplace = getappdata(hfig,'hierinplace'); + +[gIX, numU] = SqueezeGroupIX(gIX); +[C,~] = FindCentroid(gIX,M); + +str = get(hObject,'String'); +if ~isempty(str), + temp = textscan(str,'%f',1); + thres= temp{:}; + + IX_tree = clusterdata(C,'criterion','distance',... + 'distance','correlation','cutoff',thres); + + if hierinplace, + % sort to keep clusters in place + U = unique(IX_tree,'stable'); + temp = zeros(size(IX_tree)); + for i = 1:length(U), + temp(IX_tree==U(i)) = i; + end + IX_tree = temp; + end + + % update gIX + temp = zeros(size(gIX)); + for i = 1:numU, + temp(gIX==i) = IX_tree(i); + end + gIX = temp; + + cIX = getappdata(hfig,'cIX'); + UpdateIndices(hfig,cIX,gIX,length(unique(IX_tree))); + RefreshFigure(hfig); +end +end + +function checkbox_hierinplace_Callback(hObject,~) +hierinplace = get(hObject,'Value'); +hfig = getParentFigure(hObject); +setappdata(hfig,'hierinplace',hierinplace); +end + +function pushbutton_corrplot_Callback(hObject,~) +hfig = getParentFigure(hObject); +gIX = getappdata(hfig,'gIX'); +M = getappdata(hfig,'M'); +[C,~] = FindCentroid(gIX,M); +coeffs = corr(C');%corr(C(1,:)',C(2,:)') +CorrPlot(coeffs); +end + +function CorrPlot(coeffs) +im = coeffs; + +% red-white-blue colormap +cmap = zeros(64,3); +cmap(:,1) = [linspace(0,1,32), linspace(1,1,32)]; +cmap(:,2) = [linspace(0,1,32), linspace(1,0,32)]; +cmap(:,3) = [linspace(1,1,32), linspace(1,0,32)]; +minlim = -1; %min(min(im)); +maxlim = 1; %max(max(im)); + +RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap + +figure('Position',[1000,200,500,500]); +image(RGB); axis equal; axis tight; + +for i = 1:size(im,2), % horizontal.. because of image axis + for j = 1:size(im,1), + text(i-0.3, j, num2str(round(im(i,j)*100)/100));%, 'Units', 'data') + end +end +end + +function RGB = ImageToRGB(im,cmap,minlim,maxlim) +L = size(cmap,1); +ix = round(interp1(linspace(minlim,maxlim,L),1:L,im,'linear','extrap')); +RGB = reshape(cmap(ix,:),[size(ix) 3]); % Make RGB image from scaled. +end + +function pushbutton_plot4x4_Callback(hObject,~) +hfig = getParentFigure(hObject); +PlotBy16StimsFromGUI(hfig); +end + +%% ----- tab five ----- (Saved Clusters) + +%% row 1: Class + +function popup_classmenu_Callback(hObject,~) +classID = get(hObject,'Value') - 1; +if classID>0, + hfig = getParentFigure(hObject); + UpdateClassID(hfig,classID); % update popupmenu +end +end + +function edit_editclassname_Callback(hObject,~) +str = get(hObject,'String'); + +hfig = getParentFigure(hObject); +Class = getappdata(hfig,'Class'); +classID = getappdata(hfig,'classID'); +Class(classID).name = str; +setappdata(hfig,'Class',Class); +UpdateClassID(hfig,classID); % update popupmenu +end + +function pushbutton_saveclass_Callback(hObject,~) +hfig = getParentFigure(hObject); +SaveClass(hfig); +end + +function edit_class_header_Callback(hObject,~) +% get/format range +str = get(hObject,'String'); + +hfig = getParentFigure(hObject); +setappdata(hfig,'classheader',[str ':']); +end + +function edit_newclassname_Callback(hObject,~) +str = get(hObject,'String'); +hfig = getParentFigure(hObject); +setappdata(hfig,'newclassname',str); +end + +function pushbutton_makeclass_Callback(hObject,~) +hfig = getParentFigure(hObject); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); +classheader = getappdata(hfig,'classheader'); +newclassname = getappdata(hfig,'newclassname'); +SaveClass(hfig,cIX,gIX,classheader,newclassname); +end + +function pushbutton_delclass_Callback(hObject,~) +choice = questdlg('Delete current class?','','Cancel','Yes','Yes'); +if strcmp(choice,'Yes'), + hfig = getParentFigure(hObject); + Class = getappdata(hfig,'Class'); + classID = getappdata(hfig,'classID'); + disp(['delete ' num2str(classID)]); + Class(classID) = []; + setappdata(hfig,'Class',Class); + + classID = max(1,classID-1); + disp(['new classID: ' num2str(classID)]); + UpdateClassID(hfig,classID); +end +end + +%% row 2: Cluster + +function popup_clusgroupmenu_Callback(hObject,~) +hfig = getParentFigure(hObject); +clusgroupID = getappdata(hfig,'clusgroupID'); +new_clusgroupID = get(hObject,'Value'); +UpdateClusGroupID(hfig,clusgroupID,new_clusgroupID); +end + +function popup_clusmenu_Callback(hObject,~) +clusID = get(hObject,'Value') - 1; +if clusID>0, + hfig = getParentFigure(hObject); + UpdateClusID(hfig,clusID); +end +end + +function edit_editclusname_Callback(hObject,~) +str = get(hObject,'String'); + +hfig = getParentFigure(hObject); +Cluster = getappdata(hfig,'Cluster'); +clusID = getappdata(hfig,'clusID'); +Cluster(clusID).name = str; +setappdata(hfig,'Cluster',Cluster); +UpdateClusID(hfig,clusID); % update popupmenu +end + +function pushbutton_saveclus_Callback(hObject,~) +hfig = getParentFigure(hObject); +SaveCluster(hfig,'current'); +end + +function edit_clus_header_Callback(hObject,~) +% get/format range +str = get(hObject,'String'); + +hfig = getParentFigure(hObject); +setappdata(hfig,'clusheader',[str ':']); +end + +function edit_newclusname_Callback(hObject,~) +str = get(hObject,'String'); +hfig = getParentFigure(hObject); +setappdata(hfig,'newclusname',str); +end + +function pushbutton_makeclus_Callback(hObject,~) +hfig = getParentFigure(hObject); +SaveCluster(hfig,'new'); +end + +function edit_setrank_Callback(hObject,~) +str = get(hObject,'String'); +C = textscan(str,'%d'); +rank = C{:}; + +hfig = getParentFigure(hObject); +Cluster = getappdata(hfig,'Cluster'); +clusID = getappdata(hfig,'clusID'); +% insert current cluster into new position = 'rank' +if rank > 1, + temp = Cluster(clusID); + Cluster(clusID) = []; + Cluster = [Cluster(1:rank-1),temp,Cluster(rank:end)]; +else + temp = Cluster(clusID); + Cluster(clusID) = []; + Cluster = [temp,Cluster(rank:end)]; +end +setappdata(hfig,'Cluster',Cluster); +% update pointer = clusID +clusID = rank; +UpdateClusID(hfig,clusID); +end + +function edit_notes_Callback(hObject,~) +str = get(hObject,'String'); +hfig = getParentFigure(hObject); +Cluster = getappdata(hfig,'Cluster'); +clusID = getappdata(hfig,'clusID'); +Cluster(clusID).notes = str; +setappdata(hfig,'Cluster',Cluster); +end + +function pushbutton_delclus_Callback(hObject,~) +choice = questdlg('Delete current cluster?','','Cancel','Yes','Yes'); +if strcmp(choice,'Yes'), + hfig = getParentFigure(hObject); + Cluster = getappdata(hfig,'Cluster'); + clusID = getappdata(hfig,'clusID'); + disp(['delete ' num2str(clusID)]); + Cluster(clusID) = []; + setappdata(hfig,'Cluster',Cluster); + + clusID = max(1,clusID-1); + disp(['new clusID: ' num2str(clusID)]); + UpdateClusID(hfig,clusID); +end +end + +%% row 3: misc + +function edit_clusUnion_Callback(hObject,~) +disp('union processing...'); +hfig = getParentFigure(hObject); +M_0 = GetM_0(hfig); +Cluster = getappdata(hfig,'Cluster'); +str = get(hObject,'String'); +if ~isempty(str), + % get/format range + str = strrep(str,':','-'); % e.g. str= '1,3,5:8'; + C = textscan(str,'%d','delimiter',','); + m = C{:}; + range = []; + for i = 1:length(m), + if m(i)>0, + range = [range,m(i)]; %#ok + else % have '-'sign, + range = [range,m(i-1)+1:-m(i)]; %#ok + end + end + + % combine + CIX = vertcat(Cluster(range).cIX); + GIX = []; % gIX to match A + for i = 1:length(range), + GIX = [GIX; ones(length(Cluster(range(i)).gIX),1)*i]; + end + [cIX,gIX,numK] = SmartUnique(CIX,GIX,M_0(CIX,:)); + UpdateIndices(hfig,cIX,gIX,numK); + RefreshFigure(hfig); +end +disp('union complete'); +end + +function [cIX,gIX,numK] = SmartUnique(CIX,GIX,M) % input: simply concatenated groups +disp('unique based on corr coeff...'); +CTRD = FindCentroid(GIX,M); + +[uCIX,ia,~] = unique(CIX); +uGIX = GIX(ia); +if length(uCIX) == 1, + counts = length(x); +else counts = hist(CIX,uCIX); +end +IX1 = find(counts==1); +CIX1 = uCIX(IX1); % single occurence +GIX1 = uGIX(IX1); +ia = find(~ismember(CIX,CIX1)); +C = CIX(ia); % multiple copies/occurence, keep all copies +G = GIX(ia); +M_s = M(ia,:); +Coeff = zeros(length(ia),1); +for i = 1:length(ia), + reg = CTRD(G(i),:); + Coeff(i) = corr(M_s(i,:)',reg'); +end +[~,I] = sort(Coeff,'descend'); % rank by coeff, so unique ('stable') gets highest +C = C(I); % finish sorting +G = G(I); +[CIX2,ia,~] = unique(C,'stable'); % keep the order! +GIX2 = G(ia); % the chosen copy for those with multiple copies +cIX = vertcat(CIX1,CIX2); +gIX = vertcat(GIX1,GIX2); +[gIX, numK] = SqueezeGroupIX(gIX); +disp('smart union complete'); +end + +function pushbutton_newclusgroup_Callback(hObject,~) +hfig = getParentFigure(hObject); +ClusGroup = getappdata(hfig,'ClusGroup'); +ClusGroup = [ClusGroup,cell(1,1)]; +setappdata(hfig,'ClusGroup',ClusGroup); + +clusgroupID = getappdata(hfig,'clusgroupID'); +new_clusgroupID = length(ClusGroup); +UpdateClusGroupID(hfig,clusgroupID,new_clusgroupID); +end + +function pushbutton_delclusgroup_Callback(hObject,~) +choice = questdlg('Delete current Cluster-group?','','Cancel','Yes','Yes'); +if strcmp(choice,'Yes'), + hfig = getParentFigure(hObject); + ClusGroup = getappdata(hfig,'ClusGroup'); + clusgroupID = getappdata(hfig,'clusgroupID'); + ClusGroup(clusgroupID) = []; + setappdata(hfig,'ClusGroup',ClusGroup); + UpdateClusGroupID(hfig,[],max(1,clusgroupID-1)); +end +end + +%% Internal functions + +function SetDataFR(hfig,dataFR) +cIX = getappdata(hfig,'cIX'); +setappdata(hfig,'dataFR',dataFR); +if dataFR, + M_0_fluo = getappdata(hfig,'M_0_fluo'); + setappdata(hfig,'M',M_0_fluo(cIX,:)); +else % Unchecked + M_0_reg = getappdata(hfig,'M_0_reg'); + setappdata(hfig,'M',M_0_reg(cIX,:)); +end +global hdataFR; +set(hdataFR,'Value',dataFR); +end + +function SaveClass(hfig,cIX,gIX,classheader,newclassname) +if ~exist('cIX','var'), + cIX = getappdata(hfig,'cIX'); +end +if ~exist('gIX','var'), + gIX = getappdata(hfig,'gIX'); +end + +Class = getappdata(hfig,'Class'); +dataname = getappdata(hfig,'dataname'); + +if ~exist('newclassname','var'), + classID = getappdata(hfig,'classID'); +else + classID = numel(Class)+1; + Class(classID).name = [classheader newclassname]; +end + +Class(classID).cIX = cIX; +Class(classID).gIX = gIX; +Class(classID).numel = length(cIX); +U = unique(gIX); +numU = length(U); +Class(classID).numK = numU; +Class(classID).dataname = dataname; + +setappdata(hfig,'Class',Class); +UpdateClassID(hfig,classID); +disp('class saved'); +end + +function UpdateClassID(hfig,classID,norefresh) +Class = getappdata(hfig,'Class'); +numK = Class(classID).numK; +i_fish = getappdata(hfig,'i_fish'); +% save +setappdata(hfig,'classID',classID); +% update GUI +global hclassname hclassmenu; +if ishandle(hclassname), + set(hclassname,'String',Class(classID).name); + temp = MakeMenu({Class.name}); + set(hclassmenu, 'String', temp); + set(hclassmenu, 'Value', classID+1); +end + +if ~exist('norefresh','var'), + UpdateIndices(hfig,Class(classID).cIX,Class(classID).gIX,numK); + RefreshFigure(hfig); +end +global VAR; +VAR(i_fish).Class = getappdata(hfig,'Class'); +end + +function UpdateClusGroupID(hfig,clusgroupID,new_clusgroupID,norefresh) +% save/update old Cluster into ClusGroup before exiting, +% as Cluster is the variable handled in hfig but not saved elsewhere +ClusGroup = getappdata(hfig,'ClusGroup'); +Cluster = getappdata(hfig,'Cluster'); +i_fish = getappdata(hfig,'i_fish'); + +ClusGroup{clusgroupID} = Cluster; +setappdata(hfig,'ClusGroup',ClusGroup); +global VAR; +VAR(i_fish).ClusGroup = CurrentClusGroup(hfig); + +Cluster = ClusGroup{new_clusgroupID}; +setappdata(hfig,'Cluster',Cluster); +setappdata(hfig,'clusgroupID',new_clusgroupID); + +global hclusgroupmenu; +if ishandle(hclusgroupmenu), + num = length(ClusGroup); + menu = {}; for i = 1:num, menu = [menu,{num2str(i)}];end + set(hclusgroupmenu,'String', menu,'Value',new_clusgroupID); +end + +if ~exist('norefresh','var'), + if numel(Cluster) == 0, + SaveCluster(hfig,'new'); + else + clusID = 1; + UpdateClusID(hfig,clusID); + end +end +end + +function ClusGroup = CurrentClusGroup(hfig) +ClusGroup = getappdata(hfig,'ClusGroup'); +Cluster = getappdata(hfig,'Cluster'); +clusgroupID = getappdata(hfig,'clusgroupID'); +ClusGroup{clusgroupID} = Cluster; +setappdata(hfig,'ClusGroup',ClusGroup); +end + +function [Cluster,clusID] = SaveCluster(hfig,state,clusheader,name) +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); + +Cluster = getappdata(hfig,'Cluster'); +dataname = getappdata(hfig,'dataname'); + +if strcmp(state,'current'), + clusID = getappdata(hfig,'clusID'); +else %if strcmp(state,'new'), + clusID = numel(Cluster)+1; + if ~exist('name','var'), + name = getappdata(hfig,'newclusname'); + end + if ~exist('clusheader','var'), + clusheader = getappdata(hfig,'clusheader'); + end + Cluster(clusID).name = [clusheader name]; +end + +Cluster(clusID).cIX = cIX; +Cluster(clusID).gIX = gIX; +Cluster(clusID).numel = length(cIX); +U = unique(gIX); +numU = length(U); +Cluster(clusID).numK = numU; +Cluster(clusID).dataname = dataname; + +setappdata(hfig,'Cluster',Cluster); +UpdateClusID(hfig,clusID); +disp('cluster saved'); +end + +function UpdateClusID(hfig,clusID) +Cluster = getappdata(hfig,'Cluster'); + +% save +setappdata(hfig,'clusID',clusID); +% update GUI +global hclusname hclusmenu; +set(hclusname,'String',Cluster(clusID).name); +menu = MakeMenu({Cluster.name}); +set(hclusmenu,'String', menu,'Value',clusID+1); +numK = Cluster(clusID).numK; + +UpdateIndices(hfig,Cluster(clusID).cIX,Cluster(clusID).gIX,numK); +RefreshFigure(hfig); +end + +function menu = MakeMenu(name) % careful: input {Class.name}, otherwise only get 1 somehow +menu = [{'(choose)'},name]; +for j=2:length(menu),menu(j)={[num2str(j-1) ': ' menu{j}]};end +end + +function [gIX, numU] = HierClus(M,gIX,isplotfig) %#ok +[gIX, numU] = SqueezeGroupIX(gIX); +[C,~] = FindCentroid(gIX,M); +D = pdist(C,'correlation'); +tree = linkage(C,'average','correlation'); +leafOrder = optimalleaforder(tree,D); + +if numU>1, + if exist('isplotfig','var'), + figure('Position',[100 400 1300 400]); + subplot(1,3,1); + CORR = corr(C'); + imagesc(CORR);axis equal;axis tight + + subplot(1,3,2); + dendrogram(tree,numU,'orientation','right','reorder',leafOrder); + set(gca,'YDir','reverse'); + colormap(jet); + + subplot(1,3,3); + C2 = C(leafOrder,:); + CORR2 = corr(C2'); + imagesc(CORR2);axis equal;axis tight + end + % sort for uniform colorscale + temp = zeros(size(gIX)); + for i = 1:numU, + temp(gIX==leafOrder(i)) = i; % = T(i) for clusters segmented from tree + end + gIX = temp; +end +end + +function gIX = HierClusDirect(C,gIX,numU) +D = pdist(C,'correlation'); +tree = linkage(C,'average','correlation'); +leafOrder = optimalleaforder(tree,D); + +% sort for uniform colorscale +temp = zeros(size(gIX)); +for i = 1:numU, + temp(gIX==leafOrder(i)) = i; % = T(i) for clusters segmented from tree +end +gIX = temp; +end + +function [gIX, numK] = SqueezeGroupIX(gIX) +U = unique(gIX); +numK = length(U); +for i = 1:numK, + old = U(i); + gIX(gIX==old) = i; +end +end + +% frequently used, updates cell-index,group-index,cluster-number. set-operations included in here. +function UpdateIndices(hfig,cIX,gIX,numK) +global hback hopID; +M_0 = GetM_0(hfig); +if ~exist('gIX','var'), + gIX = getappdata(hfig,'gIX'); +end +if ~exist('cIX','var'), + cIX = getappdata(hfig,'cIX'); +end + +% update cache +bC = getappdata(hfig,'bCache'); +cIX_last = getappdata(hfig,'cIX'); +gIX_last = getappdata(hfig,'gIX'); +if ~(isequal(cIX_last,cIX) && isequal(gIX_last,gIX)), + bC.cIX = [cIX_last,bC.cIX]; + bC.gIX = [gIX_last,bC.gIX]; + bC.numK = [getappdata(hfig,'numK'),bC.numK]; + set(hback,'enable','on'); + if length(bC.cIX)>20, + bC.cIX(end) = []; + bC.gIX(end) = []; + bC.numK(end) = []; + end +end + +% set operations, if applicable +opID = getappdata(hfig,'opID'); +if opID ~= 0, + switch opID, + case 1, + disp('union'); + [~,ia,ib] = union(cIX_last,cIX,'stable'); + IX = vertcat(cIX_last(ia),cIX(ib));% IX; + case 2, + disp('intersect'); + [IX,ia,~] = intersect(cIX_last,cIX); + ib = []; + case 3, + disp('setdiff'); + [IX,ia] = setdiff(cIX_last,cIX); + ib = []; + case 4, + disp('rev setdiff'); + % swap sequence, then same as opID==3 + temp = cIX; + cIX = cIX_last; + cIX_last = temp; + temp = gIX; + gIX = gIX_last; + gIX_last = temp; + [IX,ia] = setdiff(cIX_last,cIX); + ib = []; + case 5, + disp('setxor'); + [IX,ia,ib] = setxor(cIX_last,cIX); + case 6, + disp('smartUnion'); + CIX = vertcat(cIX_last,cIX); + GIX = [gIX_last;gIX+max(gIX_last)]; % gIX to match + [cIX,gIX,numK] = SmartUnique(CIX,GIX,M_0(CIX,:)); + end + if opID<6, + if ~isempty(IX), + cIX = IX; + gIX = vertcat(gIX_last(ia),gIX(ib)+max(gIX_last(ia))); + numK = length(unique(gIX)); + % [gIX, numK] = SqueezeGroupIX(gIX); + else + errordlg('operation result is empty set!') + waitforbuttonpress; + end + end + set(hopID,'Value',1,'BackgroundColor',[1,1,1]); % reset + setappdata(hfig,'opID',0); +end + +if length(cIX)*size(M_0,2)<5*10^8, + % set M + M = M_0(cIX,:); + + setappdata(hfig,'M',M); + setappdata(hfig,'bCache',bC); + setappdata(hfig,'cIX',cIX); + setappdata(hfig,'gIX',gIX); + if exist('numK','var'), + setappdata(hfig,'numK',numK); + end +else + errordlg('dataset too large!') + waitforbuttonpress; +end + +% handle rankID: >=2 means write numbers as text next to colorbar +% first UpdateIndices sets rankID to 100, second sets back to 0 +rankID = getappdata(hfig,'rankID'); +if rankID>=2, + if rankID==100, + setappdata(hfig,'rankID',0); + else + setappdata(hfig,'rankID',100); + end +end +end + +% frequently used, 2 plotting functions are outside ('DrawClusters.m' and 'DrawClustersOnMap_LSh.m') +function RefreshFigure(hfig) +CInfo = getappdata(hfig,'CInfo'); +M = getappdata(hfig,'M'); +dataFR = getappdata(hfig,'dataFR'); +fictive = getappdata(hfig,'fictive'); +cIX = getappdata(hfig,'cIX'); +gIX = getappdata(hfig,'gIX'); + +numK = getappdata(hfig,'numK'); +stim = getappdata(hfig,'stim'); +anat_yx = getappdata(hfig,'anat_yx'); +anat_yz = getappdata(hfig,'anat_yz'); +anat_zx = getappdata(hfig,'anat_zx'); +isCentroid = getappdata(hfig,'isCentroid'); +clrmap = getappdata(hfig,'clrmap'); +rankscore = getappdata(hfig,'rankscore'); +rankID = getappdata(hfig,'rankID'); + +iswrite = (rankID>=2); + +if isempty(cIX), + errordlg('empty set! go back!'); + return; +end + +allAxesInFigure = findall(hfig,'type','axes'); +if ~isempty(allAxesInFigure) + delete(allAxesInFigure); +end + +figure(hfig); +% left subplot +h1 = axes('Position',[0.05, 0.04, 0.55, 0.83]); +if isCentroid, + [C,~] = FindCentroid(gIX,M); + DrawClusters(h1,C,unique(gIX),dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); +else + DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite); +end + +% right subplot +h2 = axes('Position',[0.61, 0.04, 0.35, 0.83]); +DrawClustersOnMap_LSh(CInfo,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,clrmap); +end + +function M_0 = GetM_0(hfig) +dataFR = getappdata(hfig,'dataFR'); +if dataFR, + M_0 = getappdata(hfig,'M_0_fluo'); +else + M_0 = getappdata(hfig,'M_0_reg'); +end +end + +function [C,D] = FindCentroid(gIX,M) +U = unique(gIX); +numU = length(U); +C = zeros(numU,size(M,2)); +D = zeros(numU,1); +for i=1:numU, + IX = find(gIX == U(i)); + if length(IX)==1, + C(i,:) = M(IX,:); + D(i) = 1; + else + M_s = M(IX,:); + [~,C1,~,D1] = kmeans(M_s,1,'distance','correlation'); + C(i,:) = C1; + D(i) = mean(D1); + end +end +end + +function closefigure_Callback(hfig,~) +global EXPORT_autorecover; +EXPORT_autorecover = getappdata(hfig); +end + +function fig = getParentFigure(fig) +% if the object is a figure or figure descendent, return the figure. Otherwise return []. +while ~isempty(fig) && ~strcmp('figure', get(fig,'type')) + fig = get(fig,'parent'); +end +end + +function runscript(flag_script,var_script) +switch flag_script + case 'push_cIX_gIX' + UpdateIndices(var_script{:}); + RefreshFigure(var_script{1}); +end +end +% end diff --git a/first version/GetMotorRegressor.m b/first version/GetMotorRegressor.m new file mode 100644 index 0000000..73f17e7 --- /dev/null +++ b/first version/GetMotorRegressor.m @@ -0,0 +1,86 @@ +function regressors = GetMotorRegressor(fictive) +%% generate GCaMP6f kernel +% GCaMP6f: decay half-time: 400±41; peak: 80±35 +% GCaMP6s: 1796±73, 480±24 + +fpsec = 1.97; + +tlen=size(fictive,2); +t=0:0.05:8; % sec +gc6=exp(-(t-(4+0.48))/1.8); +gc6(t<(4+0.48))=0; +t_im = 0:1/fpsec:1/fpsec*(tlen-1); +t_gc6=0:0.05:tlen; % sec +%% +% turn(data.swimStartIndT(i)-100,1) = turn_directLeft(i); %binary: left turns +% turn(data.swimStartIndT(i)-100,2) = turn_directRight(i); %binary: right turns +% turn(data.swimStartIndT(i)-100,3) = turn_direct(i); %binary: all turns (+1 for left -1 for right) +% turn(data.swimStartIndT(i)-100,4) = forward_direct(i); %binary: forward swims +% turn(data.swimStartIndT(i)-100,5) = swim_direct(i); %binary: all swims +% +% turn(data.swimStartIndT(i)-100,6) = turn_amp(i); %weighted: all turns +% turn(data.swimStartIndT(i)-100,7) = turn_left(i); %weighted: left turns +% turn(data.swimStartIndT(i)-100,8) = turn_right(i); %weighted: right turns +% turn(data.swimStartIndT(i)-100,9) = forward_amp(i); %weighted: forward swims +% turn(data.swimStartIndT(i)-100,10) = swim_amp(i); %weighted: all swims +% +% turn(data.swimStartIndT(i)-100,11) = left_amp(i); %weighted: left channel % 11,12 similar to 10?! +% turn(data.swimStartIndT(i)-100,12) = right_amp(i); %weighted: right channel +% turn(:,13) = fltCh1; %analog: left channel +% turn(:,14) = fltCh2; %analog: right channel +regressor_0={ % rows = [7,8,9,13,14]; + fictive(1,:); %weighted: right turns + fictive(2,:); %weighted: left turns + fictive(3,:); %weighted: forward swims + fictive(4,:); %analog: right channel + fictive(5,:); %analog: left channel + fictive(4,:)+fictive(5,:); %analog: average + }; +nRegType = length(regressor_0); +name_array = {'w_right','w_left','w_fwd','raw_right','raw_left','raw_all'}; + +% segment length and round off, for shuffled control +segLength = floor(tlen/80); +s4_ = tlen-mod(tlen,segLength); + +% initialize/preallocate struct +regressors(nRegType).name = []; +regressors(nRegType).im = []; +regressors(nRegType).ctrl = []; + +%% feed all regressors into struct +idx = 0; +for j=1:nRegType, %run_StimRegType_subset, + len = size(regressor_0{j},1); + for i=1:len, %run_PhotoState_subset, + idx = idx + 1; + reg = regressor_0{j}(i,:); + % idx = (j-1)*nStimRegType + i; + + regressors(idx).name = [name_array{j} '_' num2str(i)]; + % generate regressor in imaging time + regressors(idx).im = gen_reg_im(gc6, t_gc6, t_im, reg); + + % generate shuffled regressor + reg_ = reg(1:s4_); + reg2D = reshape(reg_, segLength, []); + indices = randperm(size(reg2D,2)); % reshuffle by indexing + shffreg = reg2D(:, indices); + shffreg = reshape(shffreg,1,[]); + temp = reg; + temp(1:s4_) = shffreg; + shffreg = temp; + regressors(idx).ctrl = gen_reg_im(gc6, t_gc6, t_im, shffreg); + + end +end + +end + +function reg_im = gen_reg_im(gc6, t_gc6, t_im, reg) + +temp1=interp1(t_im,reg,t_gc6,'linear','extrap'); +temp2=conv(temp1,gc6,'same'); +reg_im = interp1(t_gc6,temp2,t_im,'linear','extrap'); + +end diff --git a/GetPhotoTrans.m b/first version/GetPhotoTrans.m similarity index 100% rename from GetPhotoTrans.m rename to first version/GetPhotoTrans.m diff --git a/first version/GetStimBar.m b/first version/GetStimBar.m new file mode 100644 index 0000000..2f7215d --- /dev/null +++ b/first version/GetStimBar.m @@ -0,0 +1,25 @@ +function stimbar = GetStimBar(halfbarheight,stim) % horizontal stimulus bar. nFrames is length of actual content, barlength is raw length including padding +m1 = 0.3*ones(halfbarheight,length(stim)); +m2 = m1; % bottom half +x = stim; + +% 0 = all black; 1 = black/white; 2 = white/black; 3 = all white; 4 = all gray; +% 5 = gray/black; 6 = white/gray; 7 = black/gray; 8 = gray/white. +% 10 = forward grating (very slow, more for calibration) +% 11 = rightward grating +% 12 = leftward grating + +% top (~ projection left) +m1(:,x==0 | x==2 | x==5) = 0; % black +m1(:,x==4 | x==6 | x==7) = 0.8; % grey +m1(:,x==1 | x==3 | x==8) = 1; % white +% bottom (~ projection right) +m2(:,x==0 | x==1 | x==7) = 0; +m2(:,x==4 | x==5 | x==8) = 0.8; +m2(:,x==2 | x==3 | x==6) = 1; +% colors for OMR +stimbar = repmat(vertcat(m1,m2),[1 1 3]); % 3 color layers in dim3 +stimbar(:,x==11,1) = 1; % red ~ left +stimbar(:,x==10,2) = 1; % green ~ forward +stimbar(:,x==12,3) = 1; % blue ~ right +end diff --git a/first version/GetStimRegressor.m b/first version/GetStimRegressor.m new file mode 100644 index 0000000..ba2d9d5 --- /dev/null +++ b/first version/GetStimRegressor.m @@ -0,0 +1,186 @@ +%% +% all photoStates: +% 0 = all black; 1 = black/white; 2 = white/black; 3 = all white; 4 = all gray; +% 5 = gray/black; 6 = white/gray; 7 = black/gray; 8 = gray/white. +% 10 = forward grating (very slow, more for calibration) +% 11 = rightward grating +% 12 = leftward grating + +% 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 +% 00 01 02 03 10 11 12 13 20 21 22 23 30 31 32 33 +% transition e.g.: +% photostate: 2 3 1 3 3 0 0 1 1 0 3 2 0 2 2 1 +% phototrans: 6 11 13 7 15 12 0 1 5 4 3 14 8 2 10 9 + +%% +function [regressors, names] = GetStimRegressor(stim,fishset) + +fpsec = 1.97; % should import from data... + +%% generate GCaMP6f kernel +% GCaMP6f: decay half-time: 400±41; peak: 80±35 (ms) +halftime = 0.4; +delay = 0.08; +% GCaMP6s: 1796±73, 480±24 (ms) +% halftime = 1.8; +% delay = 0.48; + +% the 2014-15 Janelia dataset used nucleus localized GCaMP6f... + +T = 8; % kernel time length in sec + +tlen=size(stim,2); +t=0:0.05:T; % sec +gc6=exp(-(t-(T/2+delay))/halftime); +gc6(t<(T/2+delay))=0; +t_im = 0:1/fpsec:1/fpsec*(tlen-1); +t_gc6=0:0.05:tlen; % sec + +%% +if fishset == 1, + States = [0,1,2,3]; + names = {'black','phototaxis left','phototaxis right','white',... + 'right on','left on','right off','left off'}; +elseif fishset == 2, + States = [0,1,2,3,4,10,11,12]; + names = {'black','phototaxis left','phototaxis right','white','grey','OMR forward','OMR left','OMR right',... + 'left PT&OMR','right PT&OMR'}; +end +tlen=length(stim); +impulse = 6; % in frames % arbiturary at this point + +%% single photoStates +numPS = length(States); +stimPS_on = zeros(numPS,tlen); % PS = photoState +% stimPS_start = zeros(numPS,tlen); % binary, ~ impulse1 +% stimPS_stop = zeros(numPS,tlen); +for i = 1:numPS, + if ~isempty(find(stim==States(i),1)), + stimPS_on(i,:) = (stim==States(i)); + +% diff_A = [stimPS_on(i,1) diff(stimPS_on(i,:))]; +% % if i==8, diff_A(end)=-1; end % manual correction: series always ends with PS=8 +% ind_start_A = find(diff_A==1); +% ind_stop_A = find(diff_A==-1); +% +% for k = 1:length(ind_start_A), +% if (ind_start_A(k)+impulse1)=1; + + diff --git a/first version/push_cIX_gIX.m b/first version/push_cIX_gIX.m new file mode 100644 index 0000000..0663718 --- /dev/null +++ b/first version/push_cIX_gIX.m @@ -0,0 +1,10 @@ +function push_cIX_gIX(hfig,cIX,gIX,numK) +% push cIX, gIX to hfig +if ~exist('numK','var') + numK=max(gIX); +end + +var_script={hfig,cIX,gIX,numK}; +GUI_FishExplorer('_place_holder_',0,0,0,1,'push_cIX_gIX',var_script); + +end \ No newline at end of file diff --git a/first version/read_LSstack_fast_float.m b/first version/read_LSstack_fast_float.m new file mode 100644 index 0000000..b17963a --- /dev/null +++ b/first version/read_LSstack_fast_float.m @@ -0,0 +1,15 @@ +function [stack,dim]=read_LSstack_fast_float(fpath,dim) +%% reads whole stack + +info=dir(fpath); +if (isempty(info)) + error('File Does not exist'); +end + + +stack=read_LSstack_fast_float_mex64(fpath,int32(dim(1:2))); + +size(stack) + + +stack=reshape(stack,dim); diff --git a/fishpic.jpg b/fishpic.jpg new file mode 100644 index 0000000..ef40e83 Binary files /dev/null and b/fishpic.jpg differ diff --git a/read_LSstack_fast_float_mex64.mexw64 b/read_LSstack_fast_float_mex64.mexw64 new file mode 100644 index 0000000..d5df9b0 Binary files /dev/null and b/read_LSstack_fast_float_mex64.mexw64 differ