griddata3evs(V,xi,yi,zi,npiece,ln)

function [V] = griddata3evs(V,xi,yi,zi,npiece,ln)
% This function divides the volume to subvolumes and apply griddata3ev to each
% one. Each subvolume is surrounded by bigger volume that is supposed to fit when
% transformation fiels is added to grid points. The results of each
% subvolume is saved to a mat file. After datagridding of all sub volumes,
% mat files are merged to get whole volumes datagrid values.
% BUMIL: Bogazici University Medical Imaging Lab 30.09.2009

% V: Multiple volume data
% xi: scattered x coordinates
% yi: scattered y coordinates
% zi: scattered z coordinates
% npiece: # of pieces in each direction to divide the volume to subvolumes. 2 divides the volume to 8 subvolumes
% ln=5; %there are 2 sub volumes. the bigger one is bigger by ln in 3 dir.
% 25 for transformation field calculations of Alper Yaman. 

%% initialization
clc
olddir=pwd;

sz = size(xi);
% npiece = 4;
sm=round(sz/npiece);
a=sm(2);sm(2)=sm(1);sm(1)=a;sm
sm*npiece
cnt=1;
%% start datagridding process
for i=1:npiece
    for j=1:npiece
        for k=1:npiece
            %% calculate limits
            limits=[(i-1)*sm(1)+1,i*sm(1),(j-1)*sm(2)+1,j*sm(2),(k-1)*sm(3)+1,k*sm(3)];
            limitsl=[(i-1)*sm(1)+1-ln,i*sm(1)+ln,(j-1)*sm(2)+1-ln,j*sm(2)+ln,(k-1)*sm(3)+1-ln,k*sm(3)+ln];
            limitsl(limitsl<1)=1;
            if i==npiece
                limits(2)=nan;limitsl(2)=nan;
            end
            if j==npiece
                limits(4)=nan;limitsl(4)=nan;
            end
            if k==npiece
                limits(6)=nan;limitsl(6)=nan;
            end
            limits
            limitsl
            %% get sub volume
            [x,y,z,xis] = subvolume(xi,limits);
            [xl,yl,zl,xisl] = subvolume(xi,limitsl);
            yisl = subvolume(yi,limitsl);
            zisl = subvolume(zi,limitsl);
            clear vl;
            if ndims(V)==4
                for ikm=1:size(V,4)
                    vl(:,:,:,ikm)=subvolume(V(:,:,:,ikm),limitsl);
                end
            else
                vl = subvolume(V,limitsl);
            end

            %% interpolation
%             x2=xl-txl; y2=yl-tyl; z2=zl-tzl; %Actually - signs should be +, however, our T data is opposite sign.
%             t=0;clear t;xx=0;clear xx;
            %t(:,:,:,1)=xl;t(:,:,:,2)=yl;t(:,:,:,3)=zl;
            disp('interpolation started...'),tic
            v=griddata3ev(xisl,yisl,zisl,vl,x,y,z);toc
            disp(strcat('interpolation finished for ',num2str(i),'-',num2str(j),'-',num2str(k),'-',' subset...'))
            fname=strcat(num2str(cnt),'.mat');
            save(fname,'v');
            cnt=cnt+1;
            % end of interpolation
        end
    end
end
disp('Datagridding of all subvolumes are finished...');
%end datagridding process

%% merge subsets
disp('Loading subvolumes to merge');
cnt=1;
sz=size(V);
V = repmat(NaN,sz);
for i=1:npiece
    for j=1:npiece
        for k=1:npiece
            % calculate limits
            limits=[(i-1)*sm(1)+1,i*sm(1),(j-1)*sm(2)+1,j*sm(2),(k-1)*sm(3)+1,k*sm(3)];
            if i==npiece
                limits(2)=sz(2);
            end
            if j==npiece
                limits(4)=sz(1);
            end
            if k==npiece
                limits(6)=sz(3);
            end
            limits
            fname=strcat(num2str(cnt),'.mat');
            load(fname);
            if ndims(V)==3
                V(limits(3):limits(4),limits(1):limits(2),limits(5):limits(6)) = v;
            elseif ndims(V)==4
                V(limits(3):limits(4),limits(1):limits(2),limits(5):limits(6),:) = v;
            end
            cnt=cnt+1;
            delete(fname);
        end
    end
end
V(isnan(V))=0;
cd(olddir);
disp('all process finished...');