Align_LymphNode_MultiStacks.m

clc;
clear
% Aligns multiple stacks of lymph node sections, comparing 2 stacks at a time
% Sections within individual stacks are aligned but whole stacks are unaligned

% CHOOSING FILES FOR ALIGNMENT, STACK 1

waitfor (msgbox({'Select all images to be aligned within first stack.'}));

[slices1,path1] = uigetfile('*.*','Select all images to be aligned within first stack','MultiSelect','on'); % get files for first stack
image_count1 = numel(slices1); % number of images in stack to iterate through
stack1 = cell(1,numel(slices1)); % create an empty cell to contain scaled down images of stack 1
orig_rgb = cell(1,numel(slices1)); % create an empty cell to contain original images of stack 1

        % NOTE: We can improve code to not store all images -- not necessary
        % Only need to select reference imgs and maybe just load in stack
        % to be aligned

for k = 1:image_count1
    rgb_img1 = imread(fullfile(path1,slices1{k})); % read images
        % ADD BORDER
        yborder = 3000;  
        xborder = 1800; 
        %rgb_img1 = padarray(rgb_img1, [yborder, xborder], 0);
        
    orig_rgb(k) = {rgb_img1}; %storing original resolution/bordered rgb images into cell
    
    rgb_img1 = imresize(rgb_img1,0.15); % scale down images
    stack1(k) = {rgb_img1};   %storing scaled down images for alignment into cell
    
        % Current date info for file naming
        formatOut = 'mm/dd/yy';
        dte = datestr(now,formatOut);
        month = [dte(1) dte(2)];
        day = [dte(4) dte(5)];
        year = [dte(7) dte(8)];

        f_dest2 = [path1, sprintf('S2_Merge_Align_%s_%s_%s_z%d.tif',month,day,year,k)];
        %imwrite(orig_rgb{k},f_dest2);      %This just re-saves the original images
                                            %of stack 1
                                            %if you want to keep the same file naming pattern
     
end

userinput = inputdlg('Enter section number (positive integer) of desired reference image:',... % choosing fixed image for both stacks to align by
             'Ref Img for Alignment')
ref_num = str2num(userinput{1});
ref_img = stack1{ref_num};
orig_img = orig_rgb{ref_num};

% CHOOSING FILES FOR ALIGNMENT, STACK 2

waitfor (msgbox({'Select all images to be aligned within second stack.'}));

[slices2,path2] = uigetfile('*.*','Select all images to be aligned within second stack','MultiSelect','on');
image_count2 = numel(slices2);
stack2 = cell(1,numel(slices2));
orig_rgb2 = cell(1, numel(slices2)); 

for k = 1:image_count2    % loop through second stack
    rgb_img2 = imread(fullfile(path2,slices2{k}));
    % ADD BORDER
        yborder = 3000;
        xborder = 1800;  
        rgb_img2 = padarray(rgb_img2, [yborder, xborder], 0);
    
    orig_rgb2(k) = {rgb_img2}; %creating a cell of the original resolution rgb images
    
    rgb_img2 = imresize(rgb_img2,0.15);
    stack2(k) = {rgb_img2};
    %figure(k), imshow(rgb_img2);
    pause(1)
end

userinput2 = inputdlg('Enter section number (positive integer) of desired moving image:',... % choosing moving image for both stacks to align by
             'Ref Img for Alignment')
ref_num2 = str2num(userinput2{1});
ref_img2 = stack2{ref_num2};

% Turn reference/fixed image grayscale, noise filter, adjust contrast
 gfixed = rgb2gray(ref_img);
 ggfixed = medfilt2(gfixed);
 img1 = adapthisteq(ggfixed); 
 
% Turn moving image grayscale, noise filter, adjust contrast
 gmoving = rgb2gray(ref_img2);
 ggmoving = medfilt2(gmoving);
 img2 = adapthisteq(ggmoving); 

flag = true;        
while flag == true  
    
% Trace around tissue for alignment
CreateStruct.Interpreter = 'tex' ;
CreateStruct.WindowStyle = 'non-modal' ;
waitfor (msgbox ({ '\fontsize{15} In the pop-up windows that follow, a tracing tool will be used to outline the desired objects for alignment and mask out any potentially interfering artifacts.';' '; '\color{red} INSTRUCTIONS:';' ' ; '(1) Continuously hold down left mouse button and drag to trace a perimeter around object(s) of interest.';' ';'CAUTION: Releasing mouse button prematurely will terminate tracing function and could reduce accuracy of the registration.';' ';'TIP: For improved registration, trace an ROI that excludes any surrounding image artifacts, if applicable.';' ';'(2) Release mouse when finished outlining the entire object.'},'Trace Around Object(s) of Interest' , CreateStruct));
    
    f1 = figure(1); 
    imshow(img1);
    r = drawfreehand();
    mask1 = createMask(r);
    close;

    f2 = figure(2);
    imshow(img2);
    r2 = drawfreehand();
    mask2 = createMask(r2);
    close;

 % Alignment of gray reference img with chosen gray img of second stack
 
    [optimizer,metric]=imregconfig('multimodal');
     optimizer.InitialRadius = optimizer.InitialRadius/3.5;
     optimizer.MaximumIterations = 300;
    
    imrf1=imref2d(size(img1));
    imrf2=imref2d(size(img2));
    imrf3=imref2d(size(orig_img));
    
    try
        tform =  imregtform(im2double(img2).*mask2, imrf2, im2double(img1).*mask1, imrf1, 'rigid',optimizer,metric);
    catch
        warning('Error in alignment. Restarting.'); 
        %flag = false;
        continue
    end

    [optimizer,metric]=imregconfig('monomodal');
    optimizer.GradientMagnitudeTolerance = 1.00000e-04;
    optimizer.MinimumStepLength = 3.00000e-04;
    optimizer.MaximumStepLength = 6.25000e-02;
    optimizer.MaximumIterations = 300;
    optimizer.RelaxationFactor = 0.500000;
    
        try
            tform =  imregtform(im2double(img2).*mask2, imrf2, im2double(img1).*mask1, imrf1,'rigid',optimizer,metric,'InitialTransformation',tform);
                               
        catch
            warning('Error in alignment. Restarting for current slice.'); 
            %flag = false;
            continue
        end
        
    % Rescaling alignment transformation for application to original image size
        scalefactor = 6 + (2/3); %inverse of 0.15
        
        tform.T(3) = tform.T(3)*scalefactor;
        tform.T(6) = tform.T(6)*scalefactor;
        
   % Apply rescaled alignment transformation to original images
        
    for k = 1:image_count2
        [rgb_aligned3] = imwarp(orig_rgb2{k},tform,'OutputView',imrf3);

        figure(k),imshowpair(orig_img,rgb_aligned3) % Visually compare aligned second stack with fixed image of stack 1
        
        %saving files
        f_dest = [path2, sprintf('S14_Merge_Align_%s_%s_%s_z%d.tif',month,day,year,k)];
        imwrite(rgb_aligned3,f_dest,'Resolution',[300,300]); 
    end
    flag = false;
end