dynamicSPECT/matlab/extractCenters.m

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