Skip to content
This repository has been archived by the owner. It is now read-only.

Commit

Permalink
intensity values are now mean values
Browse files Browse the repository at this point in the history
MPIBR-vollrathf committed Oct 10, 2016
1 parent f0629f3 commit 1e1bd1c
Showing 4 changed files with 199 additions and 52 deletions.
90 changes: 44 additions & 46 deletions functions/CellsortApplyFilter.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean, data_g)
function cell_sig = CellsortApplyFilter(ica_segments, flims, data_g)
% cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean)
%
%CellsortApplyFilter
@@ -21,60 +21,58 @@
%

if (nargin<3)||isempty(flims)
nt = size(data_g, 3)%tiff_frames(fn);
nt = size(data_g, 3);%tiff_frames(fn);
flims = [1,nt];
else
nt = diff(flims)+1;
end
if nargin<5
subtractmean = 1;
end

[pixw,pixh] = size(data_g(:,:,1));%size(imread(fn,1));
if (nargin<4)||isempty(movm)
movm = ones(pixw,pixh);
else
movm = double(movm);
end
movm(movm==0) = 1; % Just in case there are black areas in the average image
k=0;
% k=0;

cell_sig = zeros(size(ica_segments,1), nt);
ica_segments = reshape(ica_segments, [], pixw*pixh);

fprintf('Loading %5g frames for %d ROIs.\n', nt, size(ica_segments,1))
tic
while k<nt %1:numberOfFrames
ntcurr = min(500, nt-k);
mov = zeros(pixw, pixh, ntcurr);
for j=1:ntcurr
movcurr = data_g(:,:,j+k+flims(1)-1);%imread(fn, j+k+flims(1)-1);
mov(:,:,j) = movcurr;
end
mov = (mov ./ repmat(movm, [1,1,ntcurr])) - 1; % Normalize by background and subtract mean

if subtractmean
% Subtract the mean of each frame
movtm = mean(mean(mov,1),2);
mov = mov - repmat(movtm,[pixw,pixh,1]);
end

mov = reshape(mov, pixw*pixh, ntcurr);
cell_sig(:, k+[1:ntcurr]) = ica_segments*mov;
%ica_segments = reshape(ica_segments, [], pixw*pixh);

k=k+ntcurr;
fprintf('Loaded %3.0f frames; ', k)
toc
% fprintf('Loading %5g frames for %d ROIs.\n', nt, size(ica_segments,1))
% tic
% while k<nt %1:numberOfFrames
% ntcurr = min(500, nt-k);
% mov = zeros(pixw, pixh, ntcurr);
% for j=1:ntcurr
% mov(:,:,j) = data_g(:,:,j+k+flims(1)-1);
% end
% mov = (mov ./ repmat(movm, [1,1,ntcurr])) - 1; % Normalize by background and subtract mean
%
% if subtractmean
% Subtract the mean of each frame
% movtm = mean(mean(mov,1),2);
% mov = mov - repmat(movtm,[pixw,pixh,1]);
% end
%
% mov = reshape(mov, pixw*pixh, ntcurr);
% cell_sig(:, k+[1:ntcurr]) = ica_segments*mov;
%
% k=k+ntcurr;
% fprintf('Loaded %3.0f frames; ', k)
% toc
% end
dataType = class(data_g);
data_g = reshape(data_g, pixw*pixh, size(data_g,3));
ica_segments = cast(reshape(ica_segments, size(ica_segments,1), pixw*pixh), dataType);
for seg = 1:size(ica_segments,1)
temp = bsxfun(@times, data_g, ica_segments(seg,:)');
temp2 = sum(temp,1);
quotient = size(ica_segments(seg,ica_segments(seg,:)>0),2);
temp2=temp2/quotient;
cell_sig(seg,:) = temp2;
end

function j = tiff_frames(fn)
%
% n = tiff_frames(filename)
%
% Returns the number of slices in a TIFF stack.
%
% Modified April 9, 2013 for compatibility with MATLAB 2012b
% for i = 1:size(data_g,2)
% for seg = 1:size(ica_segments,1)
% sig = data_g(:,i).*cast(squeeze(ica_segments(seg,:)),dataType)';
% temp = mean(sig(sig~=0));
% cell_sig(seg,i) = temp;
% end
% end

j = length(imfinfo(fn));
end
end
end
12 changes: 10 additions & 2 deletions functions/alternativePlot.m
Original file line number Diff line number Diff line change
@@ -26,7 +26,12 @@
'value',0, 'min',0, 'max', ceil(size(cell_sig_merged{1},2)/frameSteps),...
'SliderStep', [1 1]/(ceil(size(cell_sig_merged{1},2)/frameSteps)),'Callback', @redrawPlot,...
'units','normalized','Position',[0.6 0.01 0.3 0.03]);

imagetxt = uicontrol('units','normalized','Style','text',...
'Position',[0.1 0.04 0.3 0.03],...
'String','Showing All ROIs');
plottxt = uicontrol('units','normalized','Style','text',...
'Position',[0.6 0.04 0.3 0.03],...
'String','Showing All Frames');

function redrawPlot(~,~)
val1 = round(imageSlider.Value);
@@ -56,7 +61,8 @@ function redrawPlot(~,~)
kStart = 1;
kEnd = size(cell_sig_merged{1},2);
end

imagetxt.String = ['Showing ROI ' num2str(iStart) ' to ' num2str(iEnd) ' of ' num2str(subPlotCols)];
plottxt.String = ['Showing Frame ' num2str(kStart) ' to ' num2str(kEnd) ' of ' num2str(size(cell_sig_merged{1},2))];
drawEverything(iStart, iEnd, kStart, kEnd);
end

@@ -176,5 +182,7 @@ function drawEverything(imageStart, imageEnd, plotStart, plotEnd)
vline(frames(i)+sweepDelay+s*sweepInterval+sweepDuration, 'r');
end
end
xlabel('Frames');
ylabel('Intensity + Arbitrary Offset');
end
end
144 changes: 144 additions & 0 deletions main.asv
Original file line number Diff line number Diff line change
@@ -0,0 +1,144 @@
%% load files and do cross correlation alignment
%Choose input folder, where all trial folders are and output folder where
%shift corrected tiff stacks for each trial shall be generated. You can
%also choose a binning factor [x, y, t]. nSubIm ist the number of
%iterations the shift correction will do. The more the better, but of
%course it is more time consuming. The standard value here is 10.
inFolder = 'F:\Arbeit\testData\PCAICA';
outFolder = 'F:\Arbeit\testData\PCAICA\corr';
binningFactor = [1, 1, 2];
nSubIm = 10;
CellSortXCorr(inFolder, outFolder, nSubIm, binningFactor);
%% load files
%go to corr folder
outFolder = 'F:\Arbeit\testData\PCAICA';
%shifft correction on all frames
%[data_g, data_r, frames, fn_green, fn_red, fileNamesGreen, fileNamesRed] = CellSortLoadTrials(outFolder);
%shift correction only on trials
[data_g, data_r, frames, fn_green, fn_red, fileNamesGreen, fileNamesRed] = CellSortLoadTrialsShiftCorrOnTrials(outFolder);
%If you want to load single files, use the functions below.
% fn_green = 'data_g.tif';
% fn_red = 'data_r.tif';
% data_g = readTiff(fn_green);
% data_r = readTiff(fn_red);
%% Crop Data
%If you want to crop the data to get rid of the shift correction artifacts,
%type in the intervals of pixels in x and y directions you want to keep and
%evaluate this section. Use matVis(data_r) to look at your data.
xInterval = 40:380;
yInterval = 40:380;

data_g = data_g(xInterval,yInterval,:);
data_r = data_r(xInterval,yInterval,:);
%% 1. PCA
%calculate PCAs
nPCs = 100;
flims = [];
dSamp = [1 1]; %temporal, spatial downsampling
[mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = ...
CellsortPCA(fn_green, flims, nPCs, dSamp, [], [], data_g);
%% 2a. Choose PCs
%choose PCAs, type f/b in command line to go back or forward. type two number
%to select the range "1 ENTER 25 ENTER" would select the PCAs from 1 to 25

[PCuse] = CellsortChoosePCs(fn_green, mixedfilters, data_g);

%OR specify the PCAs to use manually:
%PCuse = [7,8,18,20,21,25,27,30,32,35,36,46];
%% 2b. Plot PC spectrum
%optionally, you can display the PC spectrum..

CellsortPlotPCspectrum(fn_green, CovEvals, PCuse, data_g)

%% 3a. ICA
%calculate ICAs, set mu to a value between 0 and 1, where 1 means pure
%temporal identification and 0 pure spatial. A value of 0.1-0.2 seems to be
%best.. make sure that the algorithm converges (see command line). If not
%run it again
nIC = length(PCuse);
mu = 0.2;

[ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, CovEvals, PCuse, mu, nIC);

%% 3b. Plot ICs
%optionally, plot ICs (you have to use at least 20.. or change the number
%on the buttom
tlims = [];
dt = 0.1;

figure(2)
CellsortICAplot('series', ica_filters, ica_sig, movm, tlims, dt, [], [], [1:size(ica_sig,1)]);

%% 4a. Segment contiguous regions within ICs
%search for indepent segments, you have to adjust these values well:
%smwidth: width of a smoothing filter (1 should be okay)
%thresh: threshold for normalized standard deviation in ICs (verey
% sensitive, standard: 2)
%arealims: limit for the detected areas in pixeks: [min max] or [min]
%plotting: 1 = yes or 0 = no.. shall a plot be shown or not?

smwidth = 1;
thresh = 1.3;
arealims = [200];
plotting = 1;

[ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting);

%% 4b. CellsortApplyFilter
%filter signal and extract regions signals
cell_sig = CellsortApplyFilter(ica_segments, flims, , data_g);

%% 5. CellsortFindspikes
%optionally, find spikes in signal
deconvtau = 0;
spike_thresh = 2;
normalization = 1;

[spmat, spt, spc] = CellsortFindspikes(ica_sig, spike_thresh, dt, deconvtau, normalization);

%% Show results
%optionally, plot signals
figure(2)
CellsortICAplot('series', ica_filters, ica_sig, movm, tlims, dt, 1, 2, 1:20, spt, spc);
%% caluclate correlation matrix for cell_sig
clusterThreshold = 0.9;
[ica_segments, cluster, subPlotCols, subPlotWidth, sortedIdx] = calcCorrEff(cell_sig, ica_segments, clusterThreshold);
%% plot
%acustic sweep parameters:
%number of frames after trial begin when first sweep occurrs
sweepDelay = 100;
%number of frames between sweeps
sweepInterval = 30;
%number of sweeps
numberOfSweeps = 5;
%number of frames a sweep lasts
sweepDuration = 10;

%number of time courses to show when slide bar is used
clusterSize = 10;
frameSteps = 20;

cell_sig_merged = alternativePlot(ica_segments, cluster, cell_sig, subPlotCols, sortedIdx, data_r, frames, fileNamesGreen, clusterSize, sweepDelay, sweepInterval, numberOfSweeps, sweepDuration, frameSteps);

%% let user deselect and merge rois
%merge -> change cluster ->recereate plots
merge = {[1,2,3], [5,6]};
delete = [8,4,11];
figure;
cluster_new = cluster;
%delete
for i = 1:size(delete,2)
cluster_new{delete(i)} = [];
end
%merge
for i = 1:size(merge, 2)
for j = 2:size(merge{i},2)
cluster_new{merge{i}(1)} = [cluster_new{merge{i}(1)}, cluster_new{merge{i}(j)}];
cluster_new{merge{i}(j)} = [];
end
end
cluster_new = cluster_new(~cellfun(@isempty, cluster_new));
subPlotCols = size(cluster_new, 2);

%redraw
cell_sig_merged = alternativePlot(ica_segments, cluster_new, cell_sig, subPlotCols, sortedIdx, data_r);
5 changes: 1 addition & 4 deletions main.m
Original file line number Diff line number Diff line change
@@ -86,10 +86,7 @@

%% 4b. CellsortApplyFilter
%filter signal and extract regions signals
subtractmean = 1;

cell_sig = CellsortApplyFilter(fn_green, ica_segments, flims, movm, subtractmean, data_g);

cell_sig = CellsortApplyFilter(ica_segments, flims, data_g);
%% 5. CellsortFindspikes
%optionally, find spikes in signal
deconvtau = 0;

0 comments on commit 1e1bd1c

Please sign in to comment.