function [order, orientation, time, recovered_page] = assemble_shreds(shred_directory)
%% load
tic
% cd(shred_directory)

strip=cell(100,1);
for i=0:99
    if i<10
        fname=['0',num2str(i)];
    else
        fname=num2str(i);
    end
    full_path=[shred_directory,'\',fname];
    strip{i+1}=im2double(imread([full_path,'.png']));
end

%% rectify and crop
theta=zeros(100,1);
ogw=zeros(100,1);
for i=1:100
    im=strip{i};
    img=rgb2gray(im);
    imb=img;
    imb(imb>0)=1; % set everything that's not 0 padding to 1
    SE=strel('square',3);
    imEdge=imdilate(imb,SE)&~imerode(imb,SE); % extract edges with morph ops
    
    %     figure; imshow(imEdge)
    
    [H,T,R]=hough(imEdge,'Theta',-1.5:0.01:1.5); % hough transform over range of possible angles
    % figure; imagesc(T,R,H)
    pk=houghpeaks(H,1);
    theta(i)=T(pk(2));
    im=imrotate(im,theta(i)); % rotate
    
    % working with rectified image from here
    img=rgb2gray(im);
    imb=img;
    imb(imb>0)=1; % set everything that's not 0 padding to 1
    % erode and dilate with vertical line
    SE=strel(ones(2001,1));
    tmp=imerode(imb,SE);
    tmp=imdilate(tmp,SE);
    ln=mean(tmp,1); % average over tall dimenison
    ln=ln>0.95; % if average is over 95%, say its part of the image
    ogw(i)=sum(ln); % width in original order
    
    h=size(im,1);
    if h<3022 % if strip is already straight up, sometimes height must be adjusted
        tmp=zeros(3022,size(im,2),3);
        tmp(1:h,:,:)=im;
        im=tmp;
    end
    
    strip{i}=im(:,ln,:); % write new image that is rectified and cropped
end


%% comparison "L2 norm" stoil
% extract left and right edges for each color channel
edges=zeros(100,size(strip{1},1),4,3); % (strip,row,L1 L2 R2 R1,RGB)
for i=1:100
    tmp=strip{i};
    edges(i,:,:,:)=tmp(:,[1:2,end-1:end],:);
end

% subtract 0.5 so RGB are centered at 0
edges=edges-0.5;


Lscoremat=zeros(100,100); % (center, left neighbor) score
Rscoremat=Lscoremat; % (center, right neighbor) score
flipLscore=Lscoremat; % (center, left neighbor flipped) score
flipRscore=Lscoremat; % (center, right neighbor flipped) score

for cent=1:100
    L=squeeze(edges(cent,:,2,:));
    R=squeeze(edges(cent,:,3,:));
    for i=1:100
        if i==cent
            Lscoremat(cent,i)=nan;
            Rscoremat(cent,i)=nan;
            flipLscore(cent,i)=nan;
            flipRscore(cent,i)=nan;
        else
            tmpL=squeeze(edges(i,:,2,:));
            tmpR=squeeze(edges(i,:,3,:));
            flipR=flip(tmpR);
            flipL=flip(tmpL);
            
            Rscoremat(cent,i)=sum(sum((R-tmpL).^2));
            Lscoremat(cent,i)=sum(sum((L-tmpR).^2));
            flipRscore(cent,i)=sum(sum((R-flipR).^2));
            flipLscore(cent,i)=sum(sum((L-flipL).^2));
        end
    end
end

%% find global min score, set as starting pair
upmat=zeros(100,100,2);
dnmat=zeros(100,100,2);

upmat(:,:,1)=Rscoremat;
upmat(:,:,2)=flipRscore;

dnmat(:,:,1)=Lscoremat;
dnmat(:,:,2)=flipLscore;

posvec=zeros(100,1);
flip_pos=zeros(100,1);
w=0*ogw;

sml=min(min(min(upmat))); % find global min
lc=find(upmat==sml);
lc=lc(1);
[a,b,c]=ind2sub(size(upmat),lc);

posvec(1)=a; w(1)=ogw(a); % save first position as left partner of global min
posvec(2)=b; w(2)=ogw(b); % save second position as right partner

flip_pos(2)=c==2;
% block out selected strip from being selected again
upmat(:,b,:)=nan;
dnmat(:,b,:)=nan;
upmat(:,a,:)=nan;
dnmat(:,a,:)=nan;

%% find best right side neighbor, iterate
for i=3:100
    prev_pos=posvec(i-1);
    prev_flip=flip_pos(i-1);
    if prev_flip==1 % if the last strip was flipped
        tmpmat=dnmat; % use the matrix corresponding to flipped entries
    else
        tmpmat=upmat; % use matrix corresponding to un-flipped entries
    end
    tmpmat=squeeze(tmpmat(prev_pos,:,:)); % only search for right neighbors too previous entry
    lc=find(tmpmat==min(min(tmpmat))); % find smallest score
    lc=lc(1);
    [b,c]=ind2sub(size(tmpmat),lc);
    posvec(i)=b;
    w(i)=ogw(b); % update width vector to reflect new order
    % update flip_pos based on if previous entry was flipped, and if this
    % entry was flipped
    if (prev_flip==0 && c==1) || (prev_flip==1 && c==2)
        flip_pos(i)=0;
    elseif (prev_flip==0 && c==2) || (prev_flip==1 && c==1)
        flip_pos(i)=1;
    end
    % block out selected strip from being selected agian
    upmat(:,b,:)=nan;
    dnmat(:,b,:)=nan;
end
%% write the new image
recon=zeros(size(edges,2),sum(w),3);
current=1;
for i=1:100
    if posvec(i)~=0 % acount for any pairs that couldn't be found (not needed anymore)
        tmp=strip{posvec(i)};
        if flip_pos(i)
            tmp=imrotate(tmp,180);
        end
        recon(:,current:current+w(i)-1,:)=tmp;
    end
    current=current+w(i);
end

% return
recovered_page=recon;
order=posvec';
orientation=flip_pos';
time=toc;

end