function fcenters=extractCenters(data,v0,sigma2,na) %extract response by segment centered around point v0, taking weights of surrounding pixels as gaussian with sigma2 %nd is number of time samples %v0 is nr by 3 matrix of centers %na is the neighborhood size (in pixels, should work great for odd na) nr=size(v0,1); fprintf("Sigma %f\n",sigma2); [~,~,~,nd]=size(data); fcenters=zeros(nr,nd); na2=int8(0.5*(na+1)); fs=zeros(nr,1); for vi=1:nr fv=v0(vi,2:4);%zero beased coordinates fi=v0(vi,1) idx0=int8(fv); idx=idx0; for i=1:na idx(1)=idx0(1)+i-na2; for j=1:na idx(2)=idx0(2)+j-na2; for k=1:na idx(3)=idx0(3)+k-na2; df=double(idx)-fv; fd=df*df'; fw=exp(-0.5*fd/sigma2); fcenters(fi,:)=fcenters(fi,:)+fw*reshape(data(idx(1)+1,idx(2)+1,idx(3)+1,:),1,[]); fs(fi)=fs(fi)+fw; end end end end qS=fs*ones(size(fcenters(1,:))); fcenters=fcenters./qS; end