Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
Latest commit 6bbd164 Mar 30, 2017 History
1 contributor

Users who have contributed to this file

720 lines (617 sloc) 31.3 KB
function [fieldPar,hand,varOut]=postSpikeFieldAnalysis(avgWF,ch,Fs,preSpikeMs,neuronNames,En,varargin)
% [fieldPar,hand,varOut]=postSpikeFieldAnalysis(avgWF,ch,Fs,preSpikeMs,neuronNames,En,varargin)
% Function purpose : Calculate distribution of post spike fields (PSF)
%
% Function recives : avgWF - average spike STAs over all electrodes in ch [Double [neurons x ch x samples]
% ch - the channel numbers of the channels in avgWFWaveform [NChannels,Time] - the raw voltage samples of all channels
% Fs - sampling frequency of the WFs
% preSpikeMs - pre spike time in avgWF
% neuronNames - names of neurons [2 x n], [channel numbers ; neuron number]
% En - electrode layout
% varargin ('property name','property value')
%
% Function give back : par - a structure of output parameters
% .classIE - I/E marker (I=2, E=3)
% lowpassWF - the sif waveforms with spikes removed
% hand - a structure of handles from generated plots
%
% Last updated : 14/12/14
%help, avgWF [neurons x ch x samples]
hand=[];fieldPar=[];
%% default variables
electrodePitch=100;
nearestNeighborsDistance=190;
postSpikeFieldStartMs=4;
postSpikeFieldEndMs=20;
preSpikePeakMs=2; %this can be larger since the exact spike time is defined by the algorithm
postSpike0CrossLimMs=20;
medianFilterLengthMs=7;
spikePeakWidthMs=1;
smartInterSmoothness=0.0001; %smoothing [0 1] - higher values fit is close to data (no low pass), 0.0000005 - more low pass
weightFunctionStdMs=7;
maxPostSpikeWidthMs=3;
stdThresholdCrossingSpikeInitiation=4;
preSpikeMinInitiationMs=1.5;
preSpikeMaxInitiationMs=0.5;
postSpikeCorrMs=10; %5
maxSIFMethod='spikeNearstNeigbohrs'; %maxLocalPeak
IEclassificationMethod='kernelProd'; %'kernelProd','delay2Crossing';
ECorrTh=[0];
ICorrTh=[0];
plotIEClass=0;
plotMaxWFAll=0;
PSFMethod='max';%'integral','max','extrapInt'
fieldPositionMethod='interpolatedMaxima';%'maxima','interpolatedMaxima','COM'
removeEdges=false;
dAngle4Plot=30;
maxFields4Plot=375;
plotFieldMapAllNeurons=false;
polarPlot=true;
plotElectrodeNames=true;
plotFieldVectors=true;
summaryPlotPerNeurons=false;
neuronIdxPolarPlot=false;
plotNeuronNumbersAllFields=false;
polarPlotRemoveOuliers=false;
polarPlotDistanceThreshold=[];
normalizeColorCode=true;
extrapolateMaxima=true;
markerSizeAllFields=15;
triangulateCellPosition = true;
cellPosition=[]; % [2 x nNeurons] correction to position based on spike shape [um]
preSpikeHPMs=2;
postSpikeHPMs=3;
classIE=true; %[true,false,vec]if false, all assumed inhibitory, can also be a vector with excitatory (2) and inhibitory (3) classifications (or 0 for require classification)
%internal variables that can be added as input
lowpassWF=[];
lowpassWFBaseline=[];
%% Output list of default variables
%print out default arguments and values if no inputs are given
if nargin==0
defaultArguments=who;
for i=1:numel(defaultArguments)
eval(['defaultArgumentValue=' defaultArguments{i} ';']);
if numel(defaultArgumentValue)==1
disp([defaultArguments{i} ' = ' num2str(defaultArgumentValue)]);
else
fprintf([defaultArguments{i} ' = ']);
disp(defaultArgumentValue);
end
end
return;
end
%% Collects all input variables
for i=1:2:length(varargin)
eval([varargin{i} '=' 'varargin{i+1};'])
end
%% Main code - general calculations
postSpikeFieldStartSamples=postSpikeFieldStartMs*Fs/1000;
postSpikeFieldEndSamples=postSpikeFieldEndMs*Fs/1000;
preSpikeSamples=preSpikeMs*Fs/1000;
preSpikePeakSamples=preSpikePeakMs*Fs/1000;
spikePeakWidthSamples=spikePeakWidthMs*Fs/1000;
maxPostSpikeWidthSamples=maxPostSpikeWidthMs*Fs/1000;
postSpikeCorrSamples=postSpikeCorrMs*Fs/1000;
weightFunctionStdSamples=weightFunctionStdMs*Fs/1000;
preSpikeMinInitiationSamples=preSpikeMinInitiationMs*Fs/1000;
preSpikeMaxInitiationSamples=preSpikeMaxInitiationMs*Fs/1000;
preSpikeHPSamples=preSpikeHPMs*Fs/1000;
postSpikeHPSamples=postSpikeHPMs*Fs/1000;
medianFilterSamples=round(medianFilterLengthMs*Fs/1000/2)*2+1; %has to be an odd number
postSpike0CrossLimSamples=postSpike0CrossLimMs*Fs/1000;
[nNeurons,nCh,nSamples]=size(avgWF);
timeVec=(1:nSamples)/Fs*1000-preSpikeMs;
%Build inverse map between electrode and location
[meshX,meshY]=meshgrid(1:size(En,1),1:size(En,2));
Xc(En(~isnan(En)))=meshX(~isnan(En))*electrodePitch;
Yc(En(~isnan(En)))=meshY(~isnan(En))*electrodePitch;
% get the channel with max spike for extimating spike remove segment
maxSpikeAmp=max( abs(avgWF(:,:,(preSpikeSamples-spikePeakWidthSamples/2):(preSpikeSamples+spikePeakWidthSamples/2))) ,[],3);
[~,pMaxSpikeElec]=max( maxSpikeAmp ,[],2);
for i=1:nCh
pNeighbors{i}=find(sqrt((Xc-Xc(i)).^2+(Yc-Yc(i)).^2)<=nearestNeighborsDistance);
end
if size(ch,1)>size(ch,2)
ch=ch';
end
%% pre-process the input waveforms
if isempty(lowpassWF) || isempty(lowpassWFBaseline)
lowpassWF=zeros(size(avgWF));
lowpassWFBaseline=zeros(size(avgWF));
pBaseline=1:(preSpikePeakSamples-preSpikeMinInitiationSamples);
preExtension=round(preSpikePeakSamples/8); %extend the detection point by a few samples
for i=1:nNeurons
fprintf('%d,',i);
%extract the initiation segment right before the spike (on nearest neigbohrs) and remove the mean of the initial part of this segment so that all segments start at 0
spikeInitiationWF=squeeze(avgWF(i,pNeighbors{pMaxSpikeElec(i)},(preSpikeSamples-preSpikePeakSamples):preSpikeSamples));%
spikeInitiationWF=bsxfun(@minus,spikeInitiationWF,mean(spikeInitiationWF(:,pBaseline),2) );
%calculate the spike onset (tr) according to where std increases rapidely over different electrode
stdProfile=std(spikeInitiationWF);
pSpikeOnset=min([preSpikePeakSamples-preSpikeMaxInitiationSamples,find(stdProfile > mean(stdProfile(pBaseline)) + stdThresholdCrossingSpikeInitiation*std(stdProfile(pBaseline)),1,'first')-1]);
pSpikeStart(i)=(preSpikeSamples-preSpikePeakSamples+pSpikeOnset-preExtension);
pSpikeSoftEnd=(preSpikeSamples+maxPostSpikeWidthSamples);
%weights for slow synaptic potential extraction
w1=ones(1,nSamples);
w1(pSpikeStart(i) : pSpikeSoftEnd)=0;
w1((pSpikeSoftEnd+1):(pSpikeSoftEnd+weightFunctionStdSamples*3))=1-exp(-( (1:weightFunctionStdSamples*3)/weightFunctionStdSamples).^2);
%weights for baseline extraction
pSpikeSoftEnd2=1200;
w2=ones(1,nSamples);
w2(pSpikeStart(i) : pSpikeSoftEnd2)=0;
w2((pSpikeSoftEnd2+1):end)=1-exp(-( (1:(nSamples-pSpikeSoftEnd2))/weightFunctionStdSamples/2).^2);
w2(1:pSpikeStart(i))=1-exp(-( (pSpikeStart(i):-1:1)/weightFunctionStdSamples/3).^2);
lowpassWF(i,:,:) = csaps(1:nSamples,squeeze(avgWF(i,:,:)),smartInterSmoothness,1:nSamples,w1);
lowpassWFBaseline(i,:,:) = csaps(1:nSamples,squeeze(lowpassWF(i,:,:)),1e-6,1:nSamples,w2);
%lowpassWFBaseline(i,:,:) = lowpassWF(i,:,:);
%lowpassWFBaseline(i,:,(pSpikeStart(i)-50):1200) = interp1([1:(pSpikeStart(i)-50) 1200:2000],squeeze(lowpassWF(i,:,[1:(pSpikeStart(i)-50) 1200:2000]))',(pSpikeStart(i)-50):1200)';
%plotting
%{
h(1)=subplot(2,3,1);
plot(timeVec,squeeze(avgWF(i,pMaxSpikeElec(i),:)));hold on;plot(timeVec,squeeze(lowpassWF(i,pMaxSpikeElec(i),:)));plot(timeVec,(w1-1)*50);plot(timeVec,(w2-1)*50);
xlabel('Time [ms]');axis tight;
spikeZoom=squeeze(avgWF(i,pNeighbors{pMaxSpikeElec(i)},(preSpikeSamples-preSpikePeakSamples):(preSpikeSamples+postSpikeFieldEndSamples)));%
spikeZoom=bsxfun(@minus,spikeZoom,mean(spikeZoom(:,200),2) );
h(2)=subplot(2,3,4);
plot(timeVec((preSpikeSamples-preSpikePeakSamples):(preSpikeSamples+postSpikeFieldEndSamples)),spikeZoom');
xlabel('Time [ms]');axis tight;
h(3)=subplot(2,3,[2 6]);
[hPlot,scaleFac]=activityTracePhysicalSpacePlot(h(3),1:120,squeeze(avgWF(i,:,:)),En,'traceColor','r','DrawElectrodeNumbers',1);hold on;
[hPlot]=activityTracePhysicalSpacePlot(h(3),1:120,squeeze(lowpassWF(i,:,:)),En,'scaleFac',scaleFac);hold on;
[hPlot]=activityTracePhysicalSpacePlot(h(3),1:120,squeeze(lowpassWFBaseline(i,:,:)),En,'scaleFac',scaleFac,'traceColor',[0.5 0.5 0.5]);
pause;
delete(h);
%}
end
end
%%
%calculate baseline substracted traces
%preBaseline=median(lowpassWF(:,:,(preSpikeSamples-preSpikePeakSamples):(preSpikeSamples-preSpikeMinInitiationSamples)),3);
%baselineSubstractedSIF=bsxfun(@minus,lowpassWF,preBaseline);
%baselineSubstractedSIF=lowpassWF-lowpassWFBaseline;
baselineSubstractedSIF=bsxfun(@minus,lowpassWF,lowpassWF(:,:,preSpikeSamples));
% get the channel with max field for classification
if strcmp(maxSIFMethod,'spikeNearstNeigbohrs')
%build extended grid
nNeighbors=1;
[nRowsTmp,nColsTmp]=size(En);
EnExt=NaN(nRowsTmp+nNeighbors*2,nColsTmp+nNeighbors*2);
EnExt(1+nNeighbors:end-nNeighbors,1+nNeighbors:end-nNeighbors)=En;
%find max amp electrode
[~,pSpikeElec]=min(avgWF(:,:,preSpikeSamples+1),[],2);
for i=1:nNeurons
[pX,pY]=find(EnExt==pSpikeElec(i));
pElecs=EnExt(pX-nNeighbors:pX+nNeighbors,pY-nNeighbors:pY+nNeighbors); %get electrodes in extended grid
pElecs=pElecs(~isnan(pElecs)); %remove NaNs
nElecs=numel(pElecs);
tmp=squeeze(baselineSubstractedSIF(i,pElecs,(preSpikeSamples+postSpikeFieldStartSamples):(preSpikeSamples+postSpikeFieldEndSamples)));
pIntersection=findfirst(tmp(:,2:end)>0 & tmp(:,1:end-1)<0, 2, 1);
pIntersection(pIntersection==0)=postSpikeFieldEndSamples-postSpikeFieldStartSamples;
sortedIntersection=sort(pIntersection);
postSpikeFieldEndSamplesNew(i)=(preSpikeSamples+postSpikeFieldStartSamples)+sortedIntersection(round(0.2*nElecs));
tmp=tmp(:,1:(postSpikeFieldEndSamplesNew(i)-(preSpikeSamples+postSpikeFieldStartSamples)));
[SIFscore]=mean( tmp ,2);
[~,pOrder]=sort(abs(SIFscore));
polarityScoreAll{i}=SIFscore(pOrder(round((nElecs*0.5):end))); %take only the high 50% of fields
polarityScore(i)=mean(polarityScoreAll{i});
polarityVote(i)=mean(sign(polarityScoreAll{i}));
pMaxField(i)=pElecs(pOrder(end));
%polarity for verification does not work well
%[polarityValidity(i)]=mean(sign( mean( abs(lowpassWF(i,pElecs,(preSpikeSamples+postSpikeFieldStartSamples):(preSpikeSamples+postSpikeFieldEndSamples)))-...
% abs(lowpassWFBaseline(i,pElecs,(preSpikeSamples+postSpikeFieldStartSamples):(preSpikeSamples+postSpikeFieldEndSamples))) ,3) ));
%{
h(1)=subplot(1,3,1:2);
[hPlot,scaleFac]=activityTracePhysicalSpacePlot(h(1),1:120,squeeze(avgWF(i,:,:)),En,'traceColor','r','DrawElectrodeNumbers',1);hold on;
[hPlot]=activityTracePhysicalSpacePlot(h(1),1:120,squeeze(lowpassWF(i,:,:)),En,'scaleFac',scaleFac);hold on;
[hPlot]=activityTracePhysicalSpacePlot(h(1),1:120,squeeze(lowpassWFBaseline(i,:,:)),En,'scaleFac',scaleFac,'traceColor',[0.5 0.5 0.5]);
h(2)=subplot(1,3,3);
plot(SIFscore);
%title(['polarity= ' num2str(polarityScore(i)), ' , Validity= ' num2str(polarityValidity(i))]);
pause;
delete(h);
%}
end
end
%% inhibitory excitatory classification
%determine which neurons to classify
%classes: 3 = excitatory, 2 = inhibitory, 1 = unclassified
if numel(classIE)==1
if classIE==0 %do not classify, but set all to be inhibitory
classIE=3*ones(1,nNeurons);
elseif classIE==1 %classify all
classIE=ones(1,nNeurons);
end %nothing happens for the case of one neuron in recording that was already clasified in the input
toClassify=(classIE==1);
else
toClassify=false(1,nNeurons);
end
pNotClassified=[];
if any(toClassify)
useScore=0;
if useScore
polarityThresh=0.5;
pExcit=find(polarityScore<-polarityThresh);
pInhib=find(polarityScore>polarityThresh);
pNotClassified=find(polarityScore>=-polarityThresh & polarityScore<=polarityThresh);
else
polarityThresh=0;
pExcit=find(polarityVote<-polarityThresh);
pInhib=find(polarityVote>polarityThresh);
pNotClassified=[];
%pNotClassified=find(polarityScore>=-polarityThresh & polarityScore<=polarityThresh);
end
fieldPar.polarityScore=polarityScore;
fieldPar.polarityVote=polarityVote;
fieldPar.classIE=ones(1,nNeurons);
fieldPar.classIE(pExcit)=3; %excitatory
fieldPar.classIE(pInhib)=2; %inhibitory
fieldPar.classIE(~toClassify)=classIE(~toClassify); %give the neurons that should not be classified their original classification
else
pExcit=find(classIE==3);
pInhib=find(classIE==2);
fieldPar.classIE=classIE;
end
if plotMaxWFAll
%define number of subplots
n=ceil(sqrt(min(maxFields4Plot,nNeurons)/3/5));%define images in a 3 x 5 ratio
xPlots=n*5;
yPlots=n*3;
nPlotPerPage=xPlots*yPlots;
cMap=lines(2);
cMap=[cMap;0 0 0;0 0 0];
f=figure('Position',[50 50 1800 900],'Visible','off');
for i=1:nNeurons
h=subaxis(f,yPlots,xPlots,i,'S',0.001,'M',0.001);
plot(timeVec,squeeze(avgWF(i,pMaxField(i),:)));hold on;
plot(timeVec,squeeze(lowpassWF(i,pMaxField(i),:)),'r');axis tight;
set(h,'XTickLabel',[],'YTick',[],'XTick',0,'TickLength',h.TickLength*5);
%text(h.XLim(2),h.YLim(2),[num2str(neuronNames(1,i)) '-' num2str(neuronNames(2,i))],'VerticalAlignment','top','HorizontalAlignment','right');
text(h.XLim(1),h.YLim(1),'*','color',cMap(4-fieldPar.classIE(i),:),'FontSize',18)
text(h.XLim(2),h.YLim(2),[num2str(i) ',' num2str(neuronNames(1,i)) '-' num2str(neuronNames(2,i))],'VerticalAlignment','top','HorizontalAlignment','right');
end
f.Visible='on';
end
%{
for i=1:nNeurons
f=figure;
h=axes;
[hPlot,scaleFac]=activityTracePhysicalSpacePlot(h,1:120,squeeze(avgWF(i,:,:)),En,'traceColor','r','DrawElectrodeNumbers',1);hold on;
[hPlot]=activityTracePhysicalSpacePlot(h,1:120,squeeze(lowpassWF(i,:,:)),En,'scaleFac',scaleFac);
title(['neuron ' num2str(neuronNames(:,i)') ', class = ' num2str(fieldPar.classIE(i))]);
pause;
delete(f);
end
%}
%% calculate post spike fields
fprintf('\nCalculating PSDs...');
pRelevantSamples=(preSpikeSamples+postSpikeFieldStartSamples):(preSpikeSamples+postSpikeFieldEndSamples);
% include the fact that each neuron has a different end time for integration
% pRelevantSamples=(preSpikeSamples+postSpikeFieldStartSamples):postSpikeFieldEndSamplesNew(i);
switch PSFMethod
case 'max'
%peak voltage normalized by pre spike peak
%fieldPar.val(pInhib,:)=max(lowpassWF(pInhib,:,pRelevantSamples),[],3)-mean(lowpassWF(pInhib,:,1:(preSpikeSamples-preSpikePeakSamples)),3);
%fieldPar.val(pExcit,:)=-(min(lowpassWF(pExcit,:,(1+preSpikeSamples):(preSpikeSamples+postSpike0CrossLimSamples)),[],3)-mean(lowpassWF(pExcit,:,1:(preSpikeSamples-preSpikePeakSamples)),3));
fieldPar.val(pInhib,:)=max(lowpassWF(pInhib,:,pRelevantSamples),[],3)-mean(lowpassWF(pInhib,:,(preSpikeSamples-spikePeakWidthSamples/2):(preSpikeSamples+spikePeakWidthSamples/2)),3);
fieldPar.val(pExcit,:)=-(min(lowpassWF(pExcit,:,(1+preSpikeSamples):(preSpikeSamples+postSpike0CrossLimSamples)),[],3)-mean(lowpassWF(pExcit,:,(preSpikeSamples-spikePeakWidthSamples/2):(preSpikeSamples+spikePeakWidthSamples/2)),3));
%set not classified the same as inhibitory
fieldPar.val(pNotClassified,:)=max(lowpassWF(pNotClassified,:,pRelevantSamples),[],3)-mean(lowpassWF(pNotClassified,:,(preSpikeSamples-spikePeakWidthSamples/2):(preSpikeSamples+spikePeakWidthSamples/2)),3);
case 'maxBaselineSubstracted'
fieldPar.val(pInhib,:)=max(baselineSubstractedSIF(pInhib,:,pRelevantSamples),[],3);
fieldPar.val(pExcit,:)=-(min(baselineSubstractedSIF(pExcit,:,(1+preSpikeSamples):(preSpikeSamples+postSpike0CrossLimSamples)),[],3));
%set not classified the same as inhibitory
fieldPar.val(pNotClassified,:)=max(baselineSubstractedSIF(pNotClassified,:,pRelevantSamples),[],3);
%{
IE=['?';'I';'E'];
pTmp=find(timeVec==0);
spikeMarker=ones(120,1)*nan(1,numel(timeVec));
spikeMarker(:,pTmp)=min(lowpassWF(:));
spikeMarker(:,pTmp+1)=max(lowpassWF(:));
for i=1:nNeurons;
h1=subplot(3,4,[1 11]);
[hPlot,scaleFac]=activityTracePhysicalSpacePlot(h1,1:120,squeeze(avgWF(i,:,:)),En,'traceColor','r');hold on;
activityTracePhysicalSpacePlot(h1,1:120,squeeze(lowpassWF(i,:,:)),En,'scaleFac',scaleFac,'DrawElectrodeNumbers',1);
activityTracePhysicalSpacePlot(h1,1:120,spikeMarker,En,'scaleFac',scaleFac,'DrawElectrodeNumbers',1,'traceColor',[0.7 0.7 0.7]);
title(['Neuron=' num2str(neuronNames(:,i)') 'index=' num2str(i) ', Max ch=' num2str(pMaxField(i)) ', C=' num2str(fieldPar.polarityVote(i))]);
h2=subplot(3,4,8);hCB=IntensityPhysicalSpacePlot(ch,fieldPar.val(i,:),En,'h',h2,'plotElectrodeNumbers',0);
title(IE(fieldPar.classIE(i)));
pause;
delete([h1 h2]);
end
%}
case 'baselineSubstractedIntegral' %!!!! Has to be rewritten to support separation between excitatory and inhibitory
%mean voltage normalized by pre spike mean
fieldPar.val(pInhib,:)=mean(baselineSubstractedSIF(pInhib,:,pRelevantSamples),3); %for inhibitory cells
fieldPar.val(pExcit,:)=-mean(baselineSubstractedSIF(pExcit,:,pRelevantSamples),3); %for inhibitory cells
fieldPar.val(pNotClassified,:)=mean(baselineSubstractedSIF(pNotClassified,:,pRelevantSamples),3);
otherwise
error('SIF calculation method not valid');
end
makeGaussianFit=0;
if makeGaussianFit
gaussFit.mX=zeros(1,nNeurons);
gaussFit.mY=zeros(1,nNeurons);
gaussFit.sX=zeros(1,nNeurons);
gaussFit.sY=zeros(1,nNeurons);
gaussFit.A=zeros(1,nNeurons);
gaussFit.Theta=zeros(1,nNeurons);
for i=1:nNeurons
[fitresult] = fmgaussfit(Xc,Yc,fieldPar.val(i,:)); %[amp, ang, sx, sy, xo, yo, zo]
gaussFit.A(i)=fitresult(1);
gaussFit.Theta(i)=fitresult(2);
gaussFit.sX(i)=fitresult(3);
gaussFit.sY(i)=fitresult(4);
gaussFit.mX(i)=fitresult(5);
gaussFit.mY(i)=fitresult(6);
end
end
if removeEdges
[~,pMax]=max(fieldPar.val,[],2);
[m,n]=size(En);
fieldPar.edgeNeurons=zeros(1,nNeurons);
for i=1:nNeurons
[pX,pY]=find(En==neuronNames(1,i));
if pX==1 || pX==n || pY==1 || pY==m
fieldPar.edgeNeurons(i)=1;
else
surroundingSquare=En(pY-1:pY+1,pX-1:pX+1);
if any(any(isnan(surroundingSquare)))
fieldPar.edgeNeurons(i)=2;
end
end
end
else
fieldPar.edgeNeurons=zeros(1,nNeurons); %set all neuron as ones not at the edge
end
fprintf('\nCalculating field peak...');
switch fieldPositionMethod
case 'interpolatedMaxima' %fits a 2D polynomial on a local grid of 9 points surrounding center
[m,n]=size(En);
Z=nan([m,n]);
%Z=zeros([m,n]);
fieldCoord=zeros(2,nNeurons);
for i=1:nNeurons
Z(sub2ind([m,n],Xc(ch)/electrodePitch,Yc(ch)/electrodePitch))=fieldPar.val(i,:);
[fieldCoord(:,i)] = peakfit2d(Z);
%p = polyFit2D(Z,XGrid,YGrid,2,2);f = polyVal2D(p,XGrid,YGrid,2,2);imagesc(f)
end
fieldPar.Xfield=fieldCoord(1,:)*electrodePitch;
fieldPar.Yfield=fieldCoord(2,:)*electrodePitch;
case 'medianCOM' %biased by array edges
%pTmp=fieldPar.val>median(fieldPar.val,2)*ones(1,nCh);
medSubstractedField=fieldPar.val-(median(fieldPar.val,2)*ones(1,nCh));
fieldPar.Xfield=(sum(bsxfun(@times,medSubstractedField,Xc),2)./sum(medSubstractedField,2))';
fieldPar.Yfield=(sum(bsxfun(@times,medSubstractedField,Yc),2)./sum(medSubstractedField,2))';
case 'maxima'
[PSF,pChPSF]=max(fieldPar.val,[],2);%location of field integral maxima
fieldPar.Xfield=Xc(ch(pChPSF));
fieldPar.Yfield=Yc(ch(pChPSF));
case 'fitGaussian'
[m,n]=size(En);
Z=nan([m,n]);
%Z=zeros([m,n]);
fieldCoord=zeros(2,nNeurons);
[YGrid,XGrid]=meshgrid(1:size(Z,1),1:size(Z,2));
for i=1:nNeurons
Z(sub2ind([m,n],Xc(ch)/electrodePitch,Yc(ch)/electrodePitch))=fieldPar.val(i,:);
[fitresult] = fmgaussfit(XGrid,YGrid,Z);
fieldCoord(:,i) = fitresult([5 6]);
end
fieldPar.Xfield=fieldCoord(1,:)*electrodePitch;
fieldPar.Yfield=fieldCoord(2,:)*electrodePitch;
case 'sumOfRegMax'
[m,n]=size(En);
%Z=min(fieldPar.val(:))*ones([m+2,n+2]);
%Z0=min(fieldPar.val(:))*ones([m,n]);
Z=zeros([m+2,n+2]);
Z0=zeros([m,n]);
fieldCoord=zeros(2,nNeurons);
[YGrid,XGrid]=meshgrid(1:size(Z,1),1:size(Z,2));
for i=1:nNeurons
%Z0(sub2ind([m,n],Xc(ch)/electrodePitch,Yc(ch)/electrodePitch))=fieldPar.val(i,:);
Z0(sub2ind([m,n],Xc(ch)/electrodePitch,Yc(ch)/electrodePitch))=fieldPar.val(i,:)-min(fieldPar.val(i,:));
Z(2:end-1,2:end-1)=Z0;
[ind] = find(imregionalmax(Z,8));
pTmp=find(Z(ind)>(fieldPar.val(i,pMaxField(i))/2));
nPeaks=numel(pTmp);
ys=zeros(nPeaks,1);xs=zeros(nPeaks,1);
for j=1:nPeaks
K=Z((XGrid(ind(pTmp(j)))-1):(XGrid(ind(pTmp(j)))+1),(YGrid(ind(pTmp(j)))-1):(YGrid(ind(pTmp(j)))+1));
% approximate polynomial parameter
a = (K(2,1)+K(1,1)-2*K(1,2)+K(1,3)-2*K(3,2)-2*K(2,2)+K(2,3)+K(3,1)+K(3,3));
b = (K(3,3)+K(1,1)-K(1,3)-K(3,1));
c = (-K(1,1)+K(1,3)-K(2,1)+K(2,3)-K(3,1)+K(3,3));
%d = (2*K(2,1)-K(1,1)+2*K(1,2)-K(1,3)+2*K(3,2)+5*K(2,2)+2*K(2,3)-K(3,1)-K(3,3));
e = (-2*K(2,1)+K(1,1)+K(1,2)+K(1,3)+K(3,2)-2*K(2,2)-2*K(2,3)+K(3,1)+K(3,3));
f = (-K(1,1)-K(1,2)-K(1,3)+K(3,1)+K(3,2)+K(3,3));
% (ys,xs) is subpixel shift of peak location relative to point (2,2)
xs(j) = (6*b*c-8*a*f)/(16*e*a-9*b^2);
ys(j) = (6*b*f-8*e*c)/(16*e*a-9*b^2);
end
fieldCoord(:,i)=[mean(XGrid(ind(pTmp))-1+xs);mean(YGrid(ind(pTmp))-1+ys)];
testPos{i}=[XGrid(ind(pTmp))-1+xs YGrid(ind(pTmp))-1+ys]'*electrodePitch;
%testPos{i}=[XGrid(ind(pTmp))-1 YGrid(ind(pTmp))-1]'*electrodePitch;
end
fieldPar.Xfield=fieldCoord(1,:)*electrodePitch;
fieldPar.Yfield=fieldCoord(2,:)*electrodePitch;
%s = regionprops(L, 'Centroid');
end
%incorporate cell position
if triangulateCellPosition && isempty(cellPosition) %run cell position estimation
avgSpkWF=avgWF(:,:, (preSpikeSamples-preSpikeHPSamples+1):(preSpikeSamples+postSpikeHPSamples) )-lowpassWF(:,:, (preSpikeSamples-preSpikeHPSamples+1):(preSpikeSamples+postSpikeHPSamples) );
[est,hest]=spikePositionEstimation(avgSpkWF,ch,preSpikeHPMs,Fs,En,fieldPar.classIE,'plot3D',0);
cellPosition(1,:)=est.X;
cellPosition(2,:)=est.Y;
%{
figure;
for i=1:nNeurons
h=axes;
[hPlot,scaleFac]=activityTracePhysicalSpacePlot(h,1:120,squeeze(avgSpkWF(i,:,:)),En,'traceColor','b','DrawElectrodeNumbers',1);hold on;title(neuronNames(:,i));
pause;
delete(h);
end
%}
elseif ~triangulateCellPosition && isempty(cellPosition) %use max spike electrode as a cell position estimator
cellPosition(1,:)=Xc(neuronNames(1,:));
cellPosition(1,:)=Yc(neuronNames(1,:));
end
%create projection vectors
X=[cellPosition(1,:);fieldPar.Xfield];
Y=[cellPosition(2,:);fieldPar.Yfield];
pTmp=isnan(fieldPar.Xfield);
X(:,pTmp)=NaN;
Y(:,pTmp)=NaN;
if size(cellPosition,2)==1 && size(cellPosition,1)~=2
error('cellPosition was not entered in the correct format');
end
mag=sqrt((X(2,:)-X(1,:)).^2 + (Y(2,:)-Y(1,:)).^2);
angle=atan2(Y(2,:)-Y(1,:),X(2,:)-X(1,:));
if ~isempty(polarPlotDistanceThreshold)
pp=find(mag>=polarPlotDistanceThreshold(1) & mag<=polarPlotDistanceThreshold(2));
mag(pp)=0;
end
pPosMagI=intersect(find(mag>0 & fieldPar.edgeNeurons==0),pInhib);
pPosMagE=intersect(find(mag>0 & fieldPar.edgeNeurons==0),pExcit);
if summaryPlotPerNeurons
IE=['?';'I';'E'];
pTmp=find(timeVec==0);
spikeMarker=ones(120,1)*nan(1,numel(timeVec));
spikeMarker(:,pTmp)=min(lowpassWF(:));
spikeMarker(:,pTmp+1)=max(lowpassWF(:));
minMaxXPos=[min(Xc) max(Xc)];
minMaxYPos=[min(Yc) max(Yc)];
f=figure('position',[10 50 1500 600]);
for i=1:nNeurons
hA(1)=subplot(2,5,[1 7]);
[hPlot,scaleFac]=activityTracePhysicalSpacePlot(hA(1),1:120,squeeze(avgWF(i,:,:)),En,'traceColor',[0.2 0.1 0.8],'gridLineWidth',0.5);hold on;
%activityTracePhysicalSpacePlot(hA(1),1:120,squeeze(lowpassWF(i,:,:)),En,'scaleFac',scaleFac,'DrawElectrodeNumbers',0,'DrawGrid',0);
%activityTracePhysicalSpacePlot(hA(1),1:120,spikeMarker,En,'scaleFac',scaleFac,'traceColor',[0.7 0.7 0.7],'gridLineWidth',0.5);
title(['Neu=' num2str(neuronNames(:,i)') ',idx=' num2str(i) ',Cls=' num2str(IE(fieldPar.classIE(i))) ',Mxch=' num2str(pMaxField(i)) ',P=' num2str(fieldPar.polarityVote(i))]);
hA(2)=subplot(2,5,[3 9]);
hA(2).Clipping='off';
F = scatteredInterpolant(Xc', Yc',fieldPar.val(i,:)');
[Xtmp,Ytmp]=meshgrid([(minMaxXPos(1)):100:(minMaxXPos(2))],[(minMaxYPos(1)):100:(minMaxYPos(2))]);
[XtmpNew,YtmpNew]=meshgrid([(minMaxXPos(1)):10:(minMaxXPos(2))],[(minMaxYPos(1)):10:(minMaxYPos(2))]);
Vq = interp2(Xtmp, Ytmp,F(Xtmp,Ytmp),XtmpNew,YtmpNew,'spline');
imagesc([(minMaxXPos(1)):10:(minMaxXPos(2))],[(minMaxYPos(1)):10:(minMaxYPos(2))],Vq);hold on;
set(hA(2),'YDir','normal');
hCB=colorbar('position',[ 0.7463 0.5943 0.0064 0.3314]);
xlabel('[\mum]');
ylabel('[\mum]');
plot(Xc,Yc,'.g')
hTmp=arrow([X(1,i);Y(1,i)]',[X(2,i);Y(2,i)],'Width',4);
plot(testPos{i}(1,:),testPos{i}(2,:),'*r');
hA(3)=subplot(2,5,5);
hCB2=IntensityPhysicalSpacePlot(ch,fieldPar.val(i,:),En,'h',hA(3),'plotElectrodeNumbers',0,'plotGridLines',0,'markerSize',50,'plotColorBar',0);hold on;
hCB2=IntensityPhysicalSpacePlot(ch,maxSpikeAmp(i,:),En,'h',hA(3),'plotElectrodeNumbers',0,'plotGridLines',0,'markerSize',25);
set(hCB2,'position',[0.9129 0.7800 0.0051 0.1457],'YTick',[]);
title('out=SIF , in=spk');
pause;
delete(hA);
end
end
%% Plotting results
if polarPlot
%prepare for plotting
f=figure('position',[100 100 500 500]);
P = panel(f);
P.pack(2,2);
P.margin=8;
angleBins=(dAngle4Plot/360/2*pi):(dAngle4Plot/360*pi):(pi*2);
if polarPlotRemoveOuliers
maximalMag=median(mag([pPosMagI pPosMagE]))+6*mad(mag([pPosMagI pPosMagE]),1);
else
maximalMag=max(mag([pPosMagI pPosMagE]));
end
%inhibitory
hand.polar(1,1)=P(1, 1).select();
hRose=rose(angle(pPosMagI),angleBins);
hRose.Color=[0.8 0.2 0.2];
XdataRose = get(hRose,'Xdata');XdataRose=reshape(XdataRose,[4,numel(XdataRose)/4]);
YdataRose = get(hRose,'Ydata');YdataRose=reshape(YdataRose,[4,numel(YdataRose)/4]);
hPatch=patch(XdataRose,YdataRose,[0.8 0.2 0.2]);
set(gca,'color','k');
%compass(U,V)
hand.polar(1,2)=P(1, 2).select();
polar(0,maximalMag,'-k');hold on; %set scale for polar plot
hTmp=polar(angle(pPosMagI),mag(pPosMagI),'.r');
if neuronIdxPolarPlot
text(hTmp.XData',hTmp.YData',num2str(neuronNames(:,pPosMagI)'),'FontSize',8);
%text(hTmp.XData',hTmp.YData',num2str(pPosMagI'),'FontSize',8);
end
%excitatory
hand.polar(2,1)=P(2, 1).select();
hRose=rose(angle(pPosMagE),angleBins);
XdataRose = get(hRose,'Xdata');XdataRose=reshape(XdataRose,[4,numel(XdataRose)/4]);
YdataRose = get(hRose,'Ydata');YdataRose=reshape(YdataRose,[4,numel(YdataRose)/4]);
hPatch=patch(XdataRose,YdataRose,[0.2 0.2 0.8]);
set(gca,'color','k');
hand.polar(2,2)=P(2, 2).select();
polar(0,maximalMag,'-k');hold on; %set scale for polar plot
hTmp=polar(angle(pPosMagE),mag(pPosMagE),'.');
if neuronIdxPolarPlot
text(hTmp.XData',hTmp.YData',num2str(neuronNames(:,pPosMagE)'),'FontSize',8);
%text(hTmp.XData',hTmp.YData',num2str(pPosMagE'),'FontSize',8);
end
end
%DSI=(prefered - (prefered+pi))/(prefered + (prefered+pi))
if plotFieldVectors
f=figure('position',[100 100 700 700]);
hand.hVec=axes;
hand.hVec.WarpToFill='off'; %to avoid error in arrow3 function
if plotElectrodeNames
hand.electrodeText=text(Xc,Yc,num2str(ch'),'fontsize',8,'Parent',hand.hVec,'horizontalAlignment','center');
xlim([min(Xc)-electrodePitch max(Xc)+electrodePitch]);
ylim([min(Yc)-electrodePitch max(Yc)+electrodePitch]);
hold(hand.hVec,'on');
end
%hQ=quiver(Xc(neuronNames(1,:)),Yc(neuronNames(1,:)),intdX,intdY,'filled','lineWidth',2,'MaxHeadSize',0.1,'color','k','MarkerSize',2,'MarkerFaceColor','k');
[tmpX,tmpY]=pol2cart(angle,50);
nInhib2Display=numel(pPosMagI);
cMapR=flipud([ones(1,60);(0:0.01:0.59);(0:0.01:0.59)]');
normColorI = ceil(min(mag(pPosMagI)./maximalMag,1).*60);
if ~isempty(pPosMagI)
hand.hArrowI=arrow3([X(1,pPosMagI);Y(1,pPosMagI)]',[X(1,pPosMagI)+tmpX(pPosMagI);Y(1,pPosMagI)+tmpY(pPosMagI)]','^r2',0.7,1);hold on;
for i=1:nInhib2Display
hand.hArrowI(i+1).FaceColor=cMapR(normColorI(i) ,:,:);
end
end
nExcit2Display=numel(pPosMagE);
cMapB=flipud([(0:0.01:0.59);(0:0.01:0.59);ones(1,60)]');
normColorE = ceil(min(mag(pPosMagE)./maximalMag,1).*60);
if ~isempty(pPosMagE)
hand.hArrowE=arrow3([X(1,pPosMagE);Y(1,pPosMagE)]',[X(1,pPosMagE)+tmpX(pPosMagE);Y(1,pPosMagE)+tmpY(pPosMagE)]','^b2',0.7,1);
for i=1:nExcit2Display
hand.hArrowE(i+1).FaceColor=cMapB(normColorE(i) ,:,:);
end
end
xlabel('X [\mum]','FontSize',14);
ylabel('Y [\mum]','FontSize',14);
end
if plotFieldMapAllNeurons
if normalizeColorCode
Ilim=0;
else
Ilim=[min(fieldPar.val(:)) max(fieldPar.val(:))];
end
n=ceil(sqrt(min(maxFields4Plot,nNeurons)/3/5));%define images in a 3 x 5 ratio
xPlots=n*5;
yPlots=n*3;
nPlotPerPage=xPlots*yPlots;
f=figure;
P = panel(f);
P.pack(yPlots,xPlots);
P.margin=0.001;
for i=1:nNeurons
hand.hAllFieldAxes(i)=P(ceil(i/xPlots),i-(ceil(i/xPlots)-1)*xPlots).select();
IntensityPhysicalSpacePlot(1:120,fieldPar.val(i,:),En,'plotElectrodeNumbers',0,'plotGridLines',0,'plotColorBar',0,'markerSize',markerSizeAllFields,'h',hand.hAllFieldAxes(i),'Ilim',Ilim);
text(Xc(neuronNames(1,i))/electrodePitch-0.5,Yc(neuronNames(1,i))/electrodePitch-0.5,'o','horizontalAlignment','center','fontsize',6);
if plotNeuronNumbersAllFields
text(0,0,num2str(i),'horizontalAlignment','left','verticalAlignment','bottom','fontsize',6);
end
line( [Xc(neuronNames(1,i)) fieldPar.Xfield(i)]/electrodePitch - 0.5 , [Yc(neuronNames(1,i)) fieldPar.Yfield(i)]/electrodePitch - 0.5 ,'color','k');
end
end
if nargout==3
varOut.lowpassWF=lowpassWF;
varOut.lowpassWFBaseline=lowpassWFBaseline;
end