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