function out = tif2bin(inPath,outPath,filePrefixName,outImageName,startSlice,numSlices,numProc)
% tif2mat
% Orion L. Kafka
% Northwestern University, Summer 2016
% A simplistic script to stack tif format grayvalue images from sectioning
% e.g. tomopy or analysis e.g. dream3d to form 3d mat files for dvc

tic
% ===== USER INPUTS =====
% Input parameters
% Set parallel processsing
%isParallel = 'true';

% Size of images
xSize = 2560; % pixels
ySize = 2560; % pixels
zSize = numSlices; % pixels (or images)

% To build file names:
startInd = startSlice;
% Ending image
endInd = startInd + zSize;
% Number of digits of postfix number
numDigits = 5;

% Output parameters
%outImageName = 'Exp021_';
outFileName = [outImageName,'.mat'];
outSmallFileName = [outImageName,'x10.mat'];

%curTime = clock;
%disp(['Start time: ',num2str(curTime(4)),'h:',num2str(curTime(5)),'m:',num2str(curTime(6)),'s::',...
%        num2str(curTime(1)),'y:',num2str(curTime(2)),'m:',num2str(curTime(3)),'d'])

% ===== START COMPUTATIONS =====

disp('Initializing variables...')
volData = zeros(xSize,ySize,zSize,'uint8'); % will actually be logical, because the output is thresholded
%countz = 1;

disp('Starting data collection...')
% ii = z, jj = x, kk = y

xSlice = 1:xSize;
ySlice = 1:ySize;

if numProc ~= 1 %strcmp(isParallel, 'true'); if there is more than one processor provided by the scheduler
    % Make a unique MATLAB parallel pool
%   pool = parpool(numProc);
    eval(['pool' outFileName ' = parpool(numProc)']);
    parfor ii = 1:zSize
        ll = ii + startInd-1;
        
	imgFile = fullfile(inPath,[filePrefixName,padWithZeros(numDigits,ll),'.tiff']);
	disp(imgFile)
	
	tifImage = imread(imgFile);

        % Read in file
        % Do some filtering and processing
        tifImage = localcontrast(tifImage,0.1,-0.1);
        tifImage = localcontrast(tifImage,0.1205,-0.205);
        tifImage = wiener2(tifImage,[4,4]);
        tifImage = medfilt2(tifImage,[5,5]);

        lvl = graythresh(tifImage);
        tifBin = imbinarize(tifImage,lvl*0.998);

        tifBin = bwareaopen(tifBin,25);
        
        
        % Save mat data
        volData(:,:,ii) = uint8(tifBin(xSlice,ySlice,:));
        
        pad = '';
        for jj = 1:numDigits
            if ll < 10^(jj)
                numPad = numDigits-jj;
                for kk = 1:numPad
                    pad = [pad,'0'];
                end
		imwrite(tifBin,...
			fullfile(outPath,[outImageName,pad,num2str(ll),'.tif']),...
			'tiff','Compression','none');
                pad = '';
                break
            end
        end
	disp(['Processing number ',num2str(ii)])
    end
    
    eval(['delete(pool' outFileName ')']); % Clean up the parallel pool
else
    for ii = 1:zSize
        ll = ii + startInd-1;
        % Build file name to read in - add proper number of padding zeros
        pad = '';
        imgFile = '';
        for jj = numDigits:-1:1
            if ll < 10^(jj)
                numPad = numDigits-jj;
                for kk = 1:numPad
                    pad = [pad,'0'];
                end
                imgFile = fullfile(inPath,[filePrefixName,pad,num2str(ll),'.tif']);
                pad = '';
            end
        end
        
        % Read in file
        tifImage = imread(imgFile,'tif');
        % Do some filtering and processing
        tifImage = localcontrast(tifImage,0.04,-1);
        lvl = graythresh(tifImage);
        tifBin = imbinarize(tifImage,lvl*1.25);

        %tifBin = bwareaopen(tifBin,50);
        
        % Save mat data
        volData(:,:,ii) = uint8(tifBin(xSlice,ySlice,:));
        % Output tif to file
        pad = '';
        for jj = 1:numDigits
            if ll < 10^(jj)
                numPad = numDigits-jj;
                for kk = 1:numPad
                    pad = [pad,'0'];
                end
                imwrite(tifBin,...
                        fullfile(outPath,[outImageName,pad,num2str(ll),'.tif']),...
                        'tiff','Compression','none');
                pad = '';
                break
            end
        end        
        disp(['Processing number ',num2str(ii),'/',num2str(zSize)])
        if ii == 5
            disp(['Estimated time = ',num2str(toc/5*zSize/60),'min']);
        end
    end
end
procTime = toc;
volSmall_10 = voxelDownsample(volData,10,'uint8');
tic
% ===== DATA OUTPUT =====
disp('Saving...')
% Save mat files
save(fullfile(outPath,outFileName),'volData','-mat','-v7.3');
save(fullfile(outPath,outSmallFileName),'volSmall_10','-mat','-v7.3');
%curTime = clock;
saveTime = toc;
%disp(['End time: ',num2str(curTime(4)),'h:',num2str(curTime(5)),'m:',num2str(curTime(6)),'s::',...
%        num2str(curTime(1)),'y:',num2str(curTime(2)),'m:',num2str(curTime(3)),'d'])
disp(['Processing time   : ',num2str((procTime)/60),'min'])
disp(['Saving time       : ',num2str((saveTime)/60),'min'])
disp(['Total elapsed time: ',num2str((procTime + saveTime)/60),'min'])
disp('Process complete.')