clear
% Enter Variables: 
z_depth = 0.300;  % Z-Stack-depth setting enter the voxel (z-step) depth
x_um = 0.57;      % x-pixel value to micron camera specific
y_um = 0.57;      % y-pixel value to micron camera specific
total_slices = 13; % z slices in each stack
total_frames = 119; % total number of time points or frames in the movie

uiwait(msgbox('Make sure your variables are accurate!'))


[filename, pathname] = uigetfile('*.tif*','File Selector');
cd(pathname);
files = dir('*.tif');
% prompt = {'Enter # of Slices:','Enter # of Frames:'};
% dlg_title = 'User Inputs';
% num_lines = 1;
% x = inputdlg(prompt,dlg_title,num_lines);
% total_slices = str2num(x{1});
% total_frames = str2num(x{2});
for kk = 1:size(files,1)
    G_file = files(kk).name;
    G_info=imfinfo(G_file);
    mImage=G_info(1).Width;
    nImage=G_info(1).Height;
    numberOfImages = length(G_info);
    G_stack=zeros(nImage,mImage,numberOfImages,'uint16');

    TifLink = Tiff(G_file, 'r');
    for i = 1:numberOfImages
        TifLink.setDirectory(i);
        G_stack(:,:,i)=TifLink.read();
    end

    % File_info = G_info.ImageDescription;
    % slice_start = regexp(File_info,'slices=','end');
    % frame_start = regexp(File_info,'frames');
    % total_slices = str2num(File_info(1,slice_start+1:frame_start-1));

    % G_time_stack = zeros(nImage,mImage,total_slices, ...
    %     numberOfImages/total_slices,'uint16');
    if numberOfImages/total_slices == total_frames
        G_time_stack = reshape(G_stack,[nImage mImage total_slices ...
            numberOfImages/total_slices]);
    else
        disp('Math Error in the Slices and Frames Values')
    end

    for iii = 1:size(G_time_stack,4)
        G_time_restack(:,:,:,iii) = double(G_time_stack(:,:,:,iii))/double(max(max(max(G_time_stack(:,:,:,iii)))));

        for ii = 1:size(G_time_stack,3)
            G_Filtered(:,:,ii,iii) = medfilt2(G_time_restack(:,:,ii,iii), [3 3]);
        end


        sub_lvl = graythresh(G_Filtered(:,:,:,iii));
        if iii==1 || isequal(fix(iii/20),iii/20) == 1
            Stack_Thresholder_v4(G_Filtered(:,:,:,iii));
            uiwait(gcf)
        end

        if F1 > 0
            GBW(:,:,:,iii) = G_Filtered(:,:,:,iii) > sub_lvl*F1;
            GBW(:,:,:,iii) = bwareaopen(GBW(:,:,:,iii),3,26);
        else
            F1 = 1.0;
            GBW(:,:,:,iii) = G_Filtered(:,:,:,iii) > sub_lvl*F1;
            GBW(:,:,:,iii) = bwareaopen(GBW(:,:,:,iii),3,26);
        end
        CC(:,:,:,iii) = bwconncomp(GBW(:,:,:,iii),26);
        Green_STATS = regionprops(CC(:,:,:,iii), G_time_stack(:,:,:,iii),{'area','MeanIntensity','Centroid','PixelValues'});

        if size(Green_STATS,1) == 1
            Green_STATS.Time_Point = iii;
            Green_STATS.Centroid_Unit = Green_STATS.Centroid;
            % pixel to unit conversion X values
            Green_STATS.Centroid_Unit(1) = Green_STATS(1, 1).Centroid(1) * x_um; 
            % pixel to unit conversion Y values
            Green_STATS.Centroid_Unit(2) = Green_STATS(1, 1).Centroid(2) * y_um; 
            % pixel to unit conversion Z values       
                if numel(Green_STATS.Centroid) < 3 
                    Green_STATS.Centroid_Unit(3) = z_depth;
                elseif numel(Green_STATS.Centroid) == 3
                    Green_STATS.Centroid_Unit(3) = ...
                        Green_STATS.Centroid(3) * z_depth; 
                end
            d = exist('Master_Stats');
            if d == 0 
                Master_Stats = Green_STATS;
            else
                Master_Stats = [Master_Stats; Green_STATS];
            end
        end
    end

    Master_ds = struct2dataset(Master_Stats);
    for i = 1:size(Master_ds,1)-1
        deltas(i,1) = pdist2(Master_ds.Centroid(i,:),Master_ds.Centroid(i+1,:));
        deltas_um(i,1) = pdist2(Master_ds.Centroid_Unit(i,:),Master_ds.Centroid_Unit(i+1,:));
        deltas(i,2) = Master_ds.Time_Point(i+1)-Master_ds.Time_Point(i);
        deltas_um(i,2) = Master_ds.Time_Point(i+1)-Master_ds.Time_Point(i);
    end
    deltas(:,3) = deltas(:,1)./deltas(:,2);
    deltas_um(:,3) = deltas_um(:,1)./deltas_um(:,2);
    Master_ds.d_Dist =  deltas(:,1);
    Master_ds.d_Time = deltas(:,2);
    Master_ds.dD_dT = deltas(:,3);
    Master_ds.d_Dist_um =  deltas_um(:,1);
    Master_ds.dD_dT_um = deltas_um(:,3);
    Master_ds.Sum_dD = sum(Master_ds.d_Dist_um);
    Master_ds.Avg_V = mean(Master_ds.dD_dT_um);
    
    figure;
    Green_Max = max(double(G_time_restack(:,:,:,1)),[],3);
    Green_Max = Green_Max/max(Green_Max(:));
    imshow(Green_Max,'InitialMagnification',800); hold on
    plot(Master_ds.Centroid(:,1),Master_ds.Centroid(:,2)); 
    title('Close figure to continue'); hold off
    uiwait(gcf)

    uisave({'Master_ds'},strrep(strrep(filename,'.tif',''),'.nd2',''))

    clear total_slices slice_start sub_lvl TifLink numberOfImages mImage ...
        nImage Master_Stats j i ii iii Green_STATS GBW CC d deltas F1 ...
        File_info frame_start G_file G_Filtered G_info G_stack ...
        G_time_restack G_time_stack pathname x prompt num_lines Green_Max
    disp('Have a GREAT day')
end