%This Script takes information from .tif files, and extracts intensity
%values from each channel of immunolabeling. It then organizes this
%information into 100 equally sized bins to determine the distribution of
%labeling within the selected region

Glomerulus=imread('A150617B5R.tif'); % reads the tiff file
Pep1=Glomerulus(:,:,1); %Determines the intensity values from channel 1;check that you have the channel you are interested in by doing figure;image(Pep1) for instance
Pep2=Glomerulus(:,:,2); %Determines the intensity values from channel 2
a=find(Pep1>254);       %The following 4 lines of code allow for matlab to read in the regions that are not squares
Pep1(a)=NaN;
b=find(Pep2>254);
Pep2(a)=NaN;


[M,N]=size(Pep1);  %Defines the variables M and N to be the number of rows (M) and the number of columns (N) for the drawn rectangle
nSub=10;           %Defines the number of bins for both rows and columns
cMat=mat2cell(Pep1,diff(round(linspace(0,M,nSub+1))),diff(round(linspace(0,N,nSub+1)))); %randomly distributes extra columns and rows when the rectangle is not equal to 10; converts from double to cell to allow different sized rows and columns

total=[];
for i=1:10
    for j=1:10
        total{i,j}=cellfun(@sum,cMat(i,j), 'UniformOutput', false); %Sums for each row;converts back to double
    end;
end;

comp=[];
for i=1:10
    for j=1:10
        for k=1:(N/10)
        comp(i,j)=sum(total{i,j}{1,:}); %Sums for each column
        end;
    end;
end;


count=[M,N];
for i=1:10
    for j=1:10
        count(i,j)=numel(cMat{i,j}); %determines the total number of pixels found for each bin
    end;
end;

zs=[];
for i=1:10
    for j=1:10
        zs{i,j}=sum((cMat{i,j}==0)); %determines the total number of pixels that are not part of picture
    end;
end;


for i=1:10
    for j=1:10
        for k=1:(N/10)
            Z3(i,j)=sum(zs{i,j}(1,:)); %Totals the number of pixels not part of picture
        end;
    end;
end;

for i=1:10
    for j=1:10
        totpix(i,j)=count(i,j)-Z3(i,j); %calculates the number of pixels used to calculate normalized intensity
    end;
end;

for i=1:10
    for j=1:10
        normintensity(i,j)=comp(i,j)/totpix(i,j); %calculates the normalized intensity
    end;
end;

l=max(normintensity);
m=max(l);
n=min(normintensity);
o=min(n);

for i=1:10
    for j=1:10
        normalizedintensity(i,j)=((normintensity(i,j)-o)/(m-o));
    end;
end;



[O,P]=size(Pep2);  %Defines the variables M and N to be the number of rows (M) and the number of columns (N) for the drawn rectangle
nSub2=10;           %Defines the number of bins for both rows and columns
cMat2=mat2cell(Pep2,diff(round(linspace(0,O,nSub+1))),diff(round(linspace(0,P,nSub+1)))); %randomly distributes extra columns and rows when the rectangle is not equal to 10; converts from double to cell to allow different sized rows and columns

total2=[];
for i=1:10
    for j=1:10
        total2{i,j}=cellfun(@sum,cMat2(i,j), 'UniformOutput', false); %Sums for each row;converts back to double
    end;
end;

comp2=[];
for i=1:10
    for j=1:10
        for k=1:(P/10)
        comp2(i,j)=sum(total2{i,j}{1,:}); %Sums for each column
        end;
    end;
end;

count2=[];
for i=1:10
    for j=1:10
        count2(i,j)=numel(cMat2{i,j}); %determines the total number of pixels found for each bin
    end;
end;

zs2=[];
for i=1:10
    for j=1:10
        zs2{i,j}=sum((cMat2{i,j}==0)); %determines the total number of pixels that are not part of picture
    end;
end;


for i=1:10
    for j=1:10
        for k=1:(P/10)
            Z32(i,j)=sum(zs2{i,j}(1,:)); %Totals the number of pixels not part of picture
        end;
    end;
end;

for i=1:10
    for j=1:10
        totpix2(i,j)=count2(i,j)-Z32(i,j); %calculates the number of pixels used to calculate normalized intensity
    end;
end;

for i=1:10
    for j=1:10
        normintensity2(i,j)=comp2(i,j)/totpix2(i,j); %calculates the normalized intensity to the total number of pixels
    end;
end;

e=max(normintensity2); %These four lines assign new variables f and h that correspond to the maximum and minimum values found in the data set
f=max(e);
g=min(normintensity2);
h=min(g);

for i=1:10
    for j=1:10
        normalizedintensity2(i,j)=((normintensity2(i,j)-h)/(f-h)); %This loop normalizes the dataset to the maximum and minimum values found in the glomerulus so that you can compare across different glomeruli
    end;
end;
        

clearvars -except normalizedintensity*;
