% This is the main program of the DVC software
%
% Developed by Soonsung Hong and Christian Franck
% California Institute of Technology, 2006-2007
%
% Updated by Jacob Notbohm
% California Institute of Technology, 2010

clear all;
clc;

% Select Incremental (inc=1) or Cumulative (inc=0) Comparison
inc = 0;

%----------------------Opening & Reading the first image-------------------
F=dir('10x_800_S1024_G575_E189_Z1_MiddleR102_0*.mat'); %d_stack
filename = F(1).name; %outputs a string with current stack file
dvc_file = load(filename);
ia = dvc_file.dec_img; % YOU'LL NEED TO UPDATE 'dec_img' TO MATCH THE VARIABLE NAME OF YOUR VOLUME STACK

%Processing Inputs:
[M,N,P]=size(ia);       % Find the size of the confocal image
w0=32;          % Subset size
d0=32;          % Subset spacing
a = 1;          % y sub-domain start
c = 1;          % x sub-domain start
b = M;       % y sub-domain end (none:M)
d = N;       % x sub-domain end (none:N)

%note: program only picks odd number of subsets
%so b=6*w0 means the height of the data set will be 5


ia = ia(a:b,c:d,:); %select sub-domain to analyze. 
%comment out above to run whole domain

%----------------------Construct 3d grid points for measurements-----------
[M,N,P]=size(ia);       % Find the size of the confocal image

x0=N/2;y0=M/2;z0=P/2;   % Find the center of the image

% Construct a set of 3d grid points centered at the center of the volume
m=round(M/2/d0)-1;n=round(N/2/d0)-1;p=round(P/2/d0)-1;
[x,y,z]=meshgrid((-n:n)*d0+x0,(-m:m)*d0+y0,(-p:p)*d0+z0);

% Initialize displacements
ua=zeros(size(x));va=ua;wa=ua;

% Preallocate storage arrays for speed
[M N P]=size(x);
uu=zeros(M,N,P,length(F)-1)*nan; vv=uu; ww=uu;

%-------------------Opening & Reading subsequent images--------------------
for stackno=1:length(F)-1
    filename = F(stackno+1).name; %outputs a string with current stack file
    dvc_file = load(filename);
    ib = dvc_file.dec_img; % YOU'LL NEED TO UPDATE 'dec_img' TO MATCH THE VARIABLE NAME OF YOUR VOLUME STACK

    %--------------------Executing DVC Code & Calculations-----------------
    % DVC translational-correlation will be run twice
    tic         % Starts computing timer (for reference only)
    % Initial guesses for translation
    ua0=ua; va0=va; wa0=wa;
    
    % Call code (run 2 iterations)
    [ub,vb,wb]=dvc_full(ia,ib,x,y,z,ua,va,wa,w0,inc);
    [uc,vc,wc]=dvc_full(ia,ib,x,y,z,ub,vb,wb,w0,inc); % 2nd iteration
    
    tm = toc;   % End time of computation of first stack. Units are sec
    disp(['Time to compute stack ',num2str(stackno),' = ',num2str(tm/3600),' hrs.'])
    disp(['Remaining Time = ',num2str(tm/3600*((length(F)-1) - stackno)),' hrs.'])
    
    % Reassign "ia" to new stack to generate incremental displacement profiles
    if inc == 1,
        ia = ib;
    end %Determines either incr. or cumulative comparison
    
    ua=uc;va=vc;wa=wc; % store current displacement fields for next increment
    
    %======================================================================
    % Update displacement results
    uu(:,:,:,stackno)=uc;
    vv(:,:,:,stackno)=vc;
    ww(:,:,:,stackno)=wc;
    % Save results after each stack
    if inc==1
        savename1=F(stackno).name; savename1=savename1(1:length(savename1)-4);
        savename2=F(stackno+1).name; savename2=savename2(1:length(savename2)-4);
        savename=['dvc_results_',savename1,'-',savename2,'.mat'];
    else
        savename1=F(1).name; savename1=savename1(1:length(savename1)-4);
        savename2=F(stackno+1).name; savename2=savename2(1:length(savename2)-4);
        savename=['dvc_results_',savename1,'-',savename2,'.mat'];
    end
    save(savename,'uc','vc','wc','w0','d0','x','y','z');
    
    disp(['Stack ',num2str(stackno),' completed.']) % Print out current stack number for reference
end
clear all;