function IntensityGraph(summary, results)

%% INTENSITY GRAPH
% Input Summary & Results excel spreadsheets from ImageJ macro to quickly
% determine the number of fusion events present in an image stack and
% identify the events as proliferating or unchanging. The function is used
% as follows:
%
% IntensityGraph('summary_sheet_name.xls', 'results_sheet_name.xls')
%
% The results are displayed via text output and the sorted data are saved
% into excel spreadsheets, e.g. "YLocations-07-Mar-2010-T12-37-47.csv"
% Figures will also display scatter plots of Area, X Location, and Y
% Location data.

warning off all

%% Open Summary & Results sheets
% Read Summary excel sheet into MATLAB
[data] = xlsread(summary);

% Store number of identified 'blobs'
num = data(:,1);
clear data

% Read Results excel sheet into MATLAB
[data] = xlsread(results);

% Store Area, X, Y data
vals = data(:,2:4);
clear data

%% Sort Area, X, Y by image number
% Create empty matrices and set initial counts
area = zeros(length(num),max(num));
X = zeros(length(num),max(num));
Y = zeros(length(num),max(num));
counta = 1;
countx = 1;
county = 1;

% For each image, find the number of blobs and sort area, X, Y data into
% rows based on the number of data points for each image. E.G. If 5 blobs
% are identified in an image, find the 5 area, X, Y data points
% corresponding to that image and place in single row - the number of which
% is the image number.
for i = 1:length(num) % number of images
    n = num(i,1); % number of blobs in each image
    d1 = vals(counta:counta+n-1,1); % pick out data points for each image
    flip1 = d1';
    area(i,:) = [flip1 zeros(1,max(num)-length(flip1))]; % store in rows
    counta = counta+n; % advance count
    d2 = vals(countx:countx+n-1,2);
    flip2 = d2';
    X(i,:) = [flip2 zeros(1,max(num)-length(flip2))];
    countx = countx+n;
    d3 = vals(county:county+n-1,3);
    flip3 = d3';
    Y(i,:) = [flip3 zeros(1,max(num)-length(flip3))];
    county = county+n;
end

%% Remove 'blob' areas under 1000 pixels
% Create new empty matrices for "large blobs" and set initial count
lgarea = zeros(size(area));
lgX = zeros(size(X));
lgY = zeros(size(Y));
count = 1;

% Any blob area under 1000 pixels isn't interesting - most cell areas are
% much larger. Go through each row and only store areas (and corresponding
% X,Y coordinates) that are larger than 1000 pixels.
for j = 1:size(X,1)
    for k = 1:size(X,2)
        if area(j,k) >= 1000
            lgarea(j,count) = area(j,k);
            lgX(j,count) = X(j,k);
            lgY(j,count) = Y(j,k);
            count = count+1;
        end
    end
    count = 1;
end

%% Remove extra columns of zeros
[r,c] = find(lgarea);
m = max(c);
newarea = lgarea(:,1:m);
newX = lgX(:,1:m);
newY = lgY(:,1:m);

%% Average nearby locations & sum nearby areas
% Create empty matrices and intial count
avgarea = zeros(size(newarea));
avgX = zeros(size(newX));
avgY = zeros(size(newY));
avgcount = 1;

% Some blobs are accidentally classified as two separate blobs, but should
% really be considered one area. If the X & Y coordinates for two areas in
% the same image are within 10 pixels of each other, the areas are added
% together and the X & Y coordinates are averaged.
for g = 1:size(newY,1)
    for h = 1:size(newY,2)-1
        if abs(newX(g,h)-newX(g,h+1))<10 && abs(newY(g,h)-newY(g,h+1))<10
            avgX(g,avgcount) = (newX(g,h)+newX(g,h+1))/2;
            avgY(g,avgcount) = (newY(g,h)+newY(g,h+1))/2;
            avgarea(g,avgcount) = newarea(g,h)+newarea(g,h+1);
            avgcount = avgcount+1;
        else
            avgX(g,avgcount) = newX(g,h);
            avgY(g,avgcount) = newY(g,h);
            avgarea(g,avgcount) = newarea(g,h);
            avgX(g,avgcount+1) = newX(g,h+1);
            avgY(g,avgcount+1) = newY(g,h+1);
            avgarea(g,avgcount+1) = newarea(g,h+1);
            avgcount = avgcount +1;
        end
    end
    avgcount = 1;
end

%% Sort data by 'blob' locations
% Create empty matrices (larger than needed)
farea = zeros(size(avgarea,1),2.*size(avgarea,2));
fX = zeros(size(avgX,1),2.*size(avgX,2));
fY = zeros(size(avgY,1),2.*size(avgY,2));

% Blobs may appear & disappear in different images. Using the Y coordinate
% data (because it is given in increasing size, generally), the coordinates
% in a column are compared to the smallest coordinate in that column. If
% they are within 150 pixels of each other, they remain in the same column
% (as they are most likely the same 'blob'). If not, all the data in their
% row is shifted to the right and a zero is placed in the column. The
% matching blob data for X and area is also moved.

for u = 1:size(fY,2)
    zY = nonzeros(avgY(:,u));
    pY = min(zY);
    
    for i = 1:size(fY,1)
        if abs(avgY(i,u)-pY)<165
            fY(i,u:size(avgY,2)) = avgY(i,u:size(avgY,2));
            fX(i,u:size(avgX,2)) = avgX(i,u:size(avgX,2));
            farea(i,u:size(avgarea,2)) = avgarea(i,u:size(avgarea,2));
        else fY(i,u) = 0; fX(i,u) = 0; farea(i,u) = 0;
            fY(i,u+1:size(avgY,2)+1) = avgY(i,u:size(avgY,2));
            fX(i,u+1:size(avgX,2)+1) = avgX(i,u:size(avgX,2));
            farea(i,u+1:size(avgarea,2)+1) = avgarea(i,u:size(avgarea,2));
        end
    end
    avgY = fY;
    avgX = fX;
    avgarea = farea;
    clear zY
end

%% Remove extra columns of zeros
[r,c] = find(fY);
m = max(c);
farea = farea(:,1:m);
fX = fX(:,1:m);
fY = fY(:,1:m);

for a = 1:size(fY,2)
    if fY(:,1)==0
        fY = fY(:,2:m);
        fX = fX(:,2:m);
        farea = farea(:,2:m);
    end
    [r,c] = find(fY);
    m = max(c);
end
clear a r c

%% Determine number of blobs/fusions & classify them
% The number of fusion events (or 'blobs') located is equal to the number
% of columns in the final area, X, and Y matrices.
s = size(farea,2);
st = num2str(s);
message = strcat('Number of fusion events=',st);
disp(message)

% Generally, fusion areas change the most drastically with proliferation.
% Thus, by finding the range of areas (max area - min area) for each fusion
% event, those with the largest ranges (>8000 pixels) are classified as
% "proliferating" and the rest are considered to be "unchanging".
lo = zeros(1,s);
hi = zeros(1,s);
for c = 1:s
    lo(1,c) = min(nonzeros(farea(:,c)));
    hi(1,c) = max(nonzeros(farea(:,c)));
end
diff = hi-lo;
pro = length(find(diff>8000));
same = s-pro;
pro = num2str(pro);
same = num2str(same);
stpro = strcat('Number of proliferating events=',pro);
stsame = strcat('Number of unchanging events=',same);
disp(stpro)
disp(stsame)

%% Plot sorted data for visual verification
% Area data for each fusion event is displayed in a scatter plot
figure;
t = 1:size(farea,1);
hold all
for c = 1:s
    plot(t,farea(:,c),'*')
end
hold off
title('\fontsize{14}Area of Fusion Events ') 
xlabel('Image Number')
ylabel('Area (pixels^2)')

% X location data for each fusion event is displayed in a scatter plot
figure;
t = 1:size(fX,1);
hold all
for c = 1:s
    plot(t,fX(:,c),'*')
end
hold off
title('\fontsize{14}X Location of Fusion Events ') 
xlabel('Image Number')
ylabel('X Location (pixels)')

% Y location data for each fusion event is displayed in a scatter plot
figure;
t = 1:size(fY,1);
hold all
for c = 1:s
    plot(t,fY(:,c),'*')
end
hold off
title('\fontsize{14}Y Location of Fusion Events ') 
xlabel('Image Number')
ylabel('Y Location (pixels)')

%% Store sorted data as new excel sheets
% For future reference, sorted data is stored into time-stamped excel
% spreadsheets.
labelY = 'YLocation-';
newlabelY = strcat(labelY, datestr(now,'dd-mmm-yyyy-THH-MM-SS'), '.xls');
xlswrite(newlabelY, fY);

labelX = 'XLocation-';
newlabelX = strcat(labelX, datestr(now,'dd-mmm-yyyy-THH-MM-SS'), '.xls');
xlswrite(newlabelX, fX);

labelarea = 'Area-';
newlabelarea = strcat(labelarea, datestr(now,'dd-mmm-yyyy-THH-MM-SS'), '.xls');
xlswrite(newlabelarea, farea);

end