[fli,corrwid,corrlen]=fimorphv2(target,n,threshtype)

function [fli,corrwid,corrlen]=fimorphv2(target,n,threshtype)

%     Thrombus Morphology
%     This technique is described in:
% 
%       T. Colace, E. Falls, X.L. Zheng, and S.L. Diamond. "Analysis of
%       morphology of platelet aggregates formed on collagen under laminar 
%       blood flow." Ann Biomed Eng. 2010.
% 
%     fimorphv2 has been used to determine the two-dimensional morphology,
%     specifically the length and width, of fluorescent platelet aggregates
%     forming on a collagen coated surface under flow. The technique uses a
%     fast, normalized two-dimensional autocorrelation procedure which is
%     performed on a stack of binary images. The method is capable of
%     performing background subtraction and selecting a threshold via a
%     stacked triangle algorithm which uses a modified version of triangle_th
%     (MATLAB file exchange: Gray image thresholding using the
%     Triangle Method by Bernard Panneton ID: 28047). The method will also
%     allow the user to provide background corrected images and specify
%     the threshold.
%
%     Inputs:
%     target - target tiff file, example: '3.tif'
%     n - number of images contained in target, example: 60
%     threshtype - 1, stacked triangle method, requires triangle_th, the user
%     will be prompted to select the correction level between -0.1 and 0.1,
%     this allows the user to correct for uneven background illumination.
%     The plot that appears shows the normalized fluorescence intensity for
%     the first, middle, and last image as a function of the correction
%     level. As the level becomes more negative the threshold value decreases
%     accepting more background. 2, the user is required to provide background
%     corrected images and is prompted to provide their ownthreshold level.
%
%     Outputs:
%     fli - an array of the fluorescence intensity of the background
%     corrected images
%     corrwid - an array of the object widths (pixels)
%     corrlen - an array of the object lengths (pixels)
%
%     plot 1 - shows the binary versions of images 1, middle, and last
%     plot 2 - shows the fluorescne intensity versus frame
%     plot 3 - shows the width/length (pixels) versus frame

delete 'temp.tif' 'temp2.tif' 'temp3.tif'
disc=50;
lehist_tot=zeros(256,1);

%Begin background correction/thresholding via triangle algorithm
if threshtype==1
    
    for i=1:n
        display(i)
        tic
        
        %Rolling ball background subtraction
        I = imread(target,'tiff',i);
        background = imopen(I,strel('disk',disc));
        I2=imsubtract(I,background);
        
        fli(i)=mean(mean(I2));
        I3=double(I2);
        %Histogram
        [lehist x]=imhist(I3/max(max(I3)));
        %Stack histogram
        lehist_tot=lehist_tot+lehist;
        %Store background subtracted images
        imwrite(I2,'temp.tif','tiff','Compression','none','Writemode','append');
        toc
    end
    
    %Triangle threshold technique
    [level,lnz]=triangle_th(lehist_tot,256,x);
    maxval=x(lnz);
    
    %User input required for illumination correction level
    for j=1:3
        img=[1 ceil(n/2) n];
        count=1;
        I2 = imread('temp.tif','tiff',img(j));
        I2= double(I2);
        for i=-.1:.01:.1
            threshval=(level+i)*maxval*max(max(I3));
            thresh=threshval/max(max(I2));
            offsetval(j,count)=i;
            
            if thresh<0
                thresh=0;
            elseif thresh>1
                thresh=1;
            end
            
            Itemp=im2bw(I2/max(max(I2)),thresh);
            Itemp = bwareaopen(Itemp,5);
            Itemp = imclose(Itemp,strel('disk',2));
            SAt(j,count)=mean(mean(Itemp));
            count=count+1;
        end
        figure(1)
        subplot(1,3,j)
        plot(offsetval(j,:),SAt(j,:))
    end
    
    i = input('Offset [-0.1 to 0.1]');
    threshval=(level+i)*maxval*max(max(I3));
    
    %Create binary images
    for i=1:n
        I2 = imread('temp.tif','tiff',i);
        I2 = double(I2);
        
        thresh=threshval/max(max(I2));
        offsetval(j,count)=i;
        
        if thresh<0
            thresh=0;
        elseif thresh>1
            thresh=1;
        end
        
        I2=im2bw(I2/max(max(I2)),thresh);
        I2 = bwareaopen(I2,5);
        I2 = imclose(I2,strel('disk',2));
        SA(i)=sum(sum(I2));
        
        imwrite(I2,'temp3.tif','tiff','Compression','none','Writemode','append');
        
        if i==1||i==ceil(n/2)||i==n
            imwrite(I2,'temp2.tif','tiff','Compression','none','Writemode','append');
        end
    end
    
%User provided background subtracted images and threshold
elseif threshtype==2
    
    thresh = input('threshold');
    
    for i=1:n
        display(i)
        tic
        
        I2 = imread(target,'tiff',i);
        I2 = double(I2);
        
        fli(i)=mean(mean(I2));
        I2=im2bw(I2/max(max(I2)),thresh/max(max(I2)));
         I2 = bwareaopen(I2,5);
         I2 = imclose(I2,strel('disk',2));
        SA(i)=sum(sum(I2));
        
        imwrite(I2,'temp.tif','tiff','Compression','none','Writemode','append');
        
        if i==1||i==ceil(n/2)||i==n
            imwrite(I2,'temp2.tif','tiff','Compression','none','Writemode','append');
        end
        
    end
    
%Threshold must be 1 or 2
else
    error('Please select threshold type 1 or type 2')
end

[x,y]=size(I2);

%Begin morphological analsysis
for i=1:n
    display(i)
    tic
    if threshtype==1
    I = imread('temp3.tif','tiff',i);
    else
    I = imread('temp.tif','tiff',i);
    end
    
    I = double(I);
    
    cenx=floor((2*x)/2);
    ceny=floor((2*y)/2);
    
    %Fast, normalized 2-d autocorrelation
    C=normxcorr2(I,I);
    
    %Blank image
    if C==0
        width=0;
        len=0;
    else
        C=C/max(max(C));
        
        %Measure length scales at 39.1% high intensity
        %Overlapping circle approximation
        D=abs(C(x,1:y)-.391);
        min(D);
        D=D-min(D);
        if isempty(find(D==0))
            width(i)=0;
        elseif length(find(D))>1
            width(i)=(y-(mean(find(D==0))))*2;
        else
            width(i)=(y-((find(D==0))))*2;
        end
        E=abs(C(1:x,y)-.391);
        E=E-min(E);
        if isempty(find(E==0))
            len(i)=0;
        elseif length(find(E))>1
            len(i)=(x-(mean(find(E==0))))*2;
        else
            len(i)=(x-((find(E==0))))*2;
        end
        %Un-normalize
        [corrwid(i)]=width(i)/(1-SA(i)/(x*y));
        [corrlen(i)]=len(i)/(1-SA(i)/(x*y));
        toc
    end
end
    
    %Plots    
    xt=[1:1:n]; %time vector
    
    figure()
    for i=1:3
        subplot(1,3,i)
        imshow(imread('temp2.tif',i))
    end
    
    figure()
    h=plot(xt,fli,'.');
    xlabel('Frame','FontSize',45,'FontWeight','bold');
    ylabel('FI','FontSize',45,'FontWeight','bold');
    set(gca,'TickLength',[.02,.025],'FontSize',40,'LineWidth',5.0,'FontWeight','bold');
    set(gca,'XLim',[0 n],'YLim',[0 300]);
    set(h,'LineStyle','none','MarkerSize',60,'Color',[0 0 0])
    
    
    figure()
    h=plot(xt,corrwid,'.','color','k','MarkerSize',60);
    hold on
    g=plot(xt,corrlen,'s','color','k','LineWidth',5,'MarkerSize',20);
    xlabel('Frame','FontSize',45,'FontWeight','bold');
    ylabel('pixel','FontSize',45,'FontWeight','bold');
    set(gca,'TickLength',[.02,.025],'FontSize',40,'LineWidth',5.0,'FontWeight','bold');
    set(gca,'XLim',[0 n],'YLim',[0 200]);
    l1 = legend([h,g],'Width','Length','Location','NorthWest','visible','off');
    legend boxoff;