function out = rotateAndCrop(inPath,outPath,runNum,numProc,varargin)
% tif2mat
% Orion L. Kafka, adapted from a script by Mouad
% Northwestern University, Winter 2018
% Crop and rotate xctImage

tic
% ===== USER INPUTS =====

imgFiles=dir(fullfile(inPath,'*.tif'));
%disp(length(imgFiles));

expNumStr = padWithZeros(3,runNum);
expNumStr = ['Exp',expNumStr];

if nargin == 6
    centroidPath = varargin{1}
    imOutPath = varargin{2}
end

% Size of images
xSize = 2560; % pixels
ySize = 2560; % pixels
% Process all the binary images
zSize = length(imgFiles); % pixels (or images)
%zSize = 150;

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

% Find the rotation angle from one of the images (it only has to be "close
% enough").
tifImage = imread(fullfile(inPath,imgFiles(floor(length(imgFiles)/6)).name));
p = regionprops(tifImage,'Extrema','Area');
[~,idx] = max([p.Area]); % Get the biggest area
sides = p(idx).Extrema(4,:) - p(idx).Extrema(6,:); % Returns dif of Bottom-Right and Left-Bottom
disp(['Sides = ',num2str(sides)]);
rotAngle = rad2deg(atan(-sides(2)/sides(1))); % The 'minus' sign compensates for the inverted y-values in image coordinates
if abs(sides(2)) > abs(sides(1))
    rotAngle = rotAngle-90;
end
disp(['Rotation angle: ',num2str(rotAngle)])

% "orientation" uses a fitted ellipsoid -- not reliable for near-square-shaped things
%rotAngle = regionprops(tifImage,'Orientation');
%rotAngle = rotAngle.Orientation

if numProc ~= 1
    % Make a unique MATLAB parallel pool
    %eval(['pool' num2str(runNum) ' = parpool(numProc);']);
    %volData = uint8(zeros(xSize,ySize,zSize));
    parfor ii = 1:zSize
%       disp(fullfile(inPath,imgFiles(ii).name))
        tifImage = imread(fullfile(inPath,imgFiles(ii).name),'tif');
        %volData(:,:,ii) = tifImage;
        tifImage = imrotate(tifImage,-rotAngle);

        [rows, columns] = find(tifImage);
        row1 = min(rows);
        row2 = max(rows);
        col1 = min(columns);
        col2 = max(columns);
        tifImage = tifImage(row1:row2, col1:col2);

        imadata2 = imfill(tifImage,'holes');
        [C H W] = FindLargestRectangles(imadata2, [1 5 1]);
        [tmp pos] = max(C(:));
        [r c] = ind2sub(size(C), pos);
        rect{ii} = [c+col1,r+row1,W(r,c),H(r,c)];
    end
    rectSmallest = [1,1,10000,10000];
    for ii = 1:zSize-1
        if rect{ii}(1) > rectSmallest(1)
            rectSmallest(1) = rect{ii}(1);
        end
        if rect{ii}(2) > rectSmallest(2)
            rectSmallest(2) = rect{ii}(2);
        end
        if rect{ii}(3) < rectSmallest(3)
            rectSmallest(3) = rect{ii}(3);
        end
        if rect{ii}(4) < rectSmallest(4)
            rectSmallest(4) = rect{ii}(4);
        end
    end
    disp(['Cropping after rotation with: xmin = ',num2str(rectSmallest(1)),...
                                      ', ymin = ',num2str(rectSmallest(2)),...
                                     ', width = ',num2str(rectSmallest(3)),...
                                    ', height = ',num2str(rectSmallest(4))]);
%    volData = zeros(rectSmallest(4)+1,rectSmallest(3)+1,zSize);
%    disp(['zsize = ',num2str(zSize)])
%    disp(['vol data size = ',num2str(size(volData))]);
    if nargin == 6
        volData = ringFill(volData,imgFiles,centroidPath,imOutPath,expNumStr,numProc);
    end

    parfor ii = 1:zSize
        tifImage = imread(fullfile(inPath,imgFiles(ii).name),'tif');
        %tifImage = volData(:,:,ii)%imread(fullfile(inPath,imgFiles(ii).name));
        tifImage = imrotate(tifImage,-rotAngle);
        tifImage = imcrop(tifImage, rectSmallest);
%        disp(['tif img size = ',num2str(size(tifImage))])
        imwrite(tifImage,fullfile(outPath,imgFiles(ii).name),'tiff','Compression','none');
    end
    save(fullfile(outPath,expNumStr),'volData','-mat','-v7.3');
    
%    if nargin == 6
%        volData = ringFill(volData,imgFiles,centroidPath,imOutPath,expNumStr,numProc);
%        save(fullfile(outPath,[expNumStr,'_fill']),'volData','-mat','-v7.3');
%    end
    eval(['delete(pool' num2str(runNum) ')']); % Clean up the parallel pool
else
    disp('Serial version not implemented')
end
procTime = toc;

disp(['Processing time   : ',num2str((procTime)/60),'min'])
disp('Process complete.')