Matlab Code: analyze collagen structure in renal cell carcinoma (RCC) (Chromophore vs Oncocytoma)

% --------------------------------- analyzeCollagen.m ------------------------------------%
%% initialization

ch_folder = {'S12-19756C8','S12-25601A6','S13-6151A4','S13-10481C8','S13-29859A9','S13-29883A8','S13-33982B4','S13-36836C5','S14-7163A3','S14-9425B4','S14-17795B2'};
on_folder = {'S12-20431B4','S12-28537A2','S12-38438B6','S13-10677A5','S13-12911A2','S13-19695B2','S13-19695B3','S13-37948A2','S14-12554B2','S14-19568A3','S14-20995C4'};
mix_folder = {'S13-3375B1'};
    
ch_fileNumber = [4 4 6 4 4 4 4 3 4 5 5];
on_fileNumber = [4 4 4 4 5 4 4 4 4 4 4];
mix_fileNumber = [4];

ch_fileCode = {
'Ta','Tb','Tc','Td', ...
'Ta','Tb','Tc','Td', ...
'T1a','T1b','T2a','T2b','T2c','T2d', ...
'Ta','Tb','Tc','Td', ...
'T1a','T1b','T1c','T1d', ...
'Ta','Tb','Tc','Td', ...
'T1b','T1c','T1d','T2d', ...
'Tb','Tc','Td', ...
'Ta','Tb','Tb2','Tc', ...
'T1a','T1b','T1c','T1d','T2a', ... %S14-9425B4
'T1a','T1b','T2a','T2b','T2d'
};

on_fileCode = {
'Ta','Ta2','Tb','Tc', ...
'Ta','Ta2','Ta3','Ta4', ...
'Ta','Ta2','Tb','Tc', ... 
'Ta','Tb','Tc','Td', ...
'T1a','T1b','T1c','T2a','T2b', ...
'T1a','T1b','T2a','T2b', ...
'T1a','T1b','T2a','T2b', ...
'T1a','T1b','T2a','T2b', ...
'T1a','T1b','T2a','T2b', ...
'Ta','Tb','Tc','Td', ...
'T1a','T1b','T2a','T2b'
};

mix_fileCode = {'T1a','T1b','T1c','T1d'};

for i = 1:length(ch_fileCode)    
    ch_fileName{i} = ['25x_',ch_fileCode{i},'_power2_stepsize1_PMT700_zoom2_kalman3_Ch1.stk'];
    ch_fileName1{i} = ['25x_',ch_fileCode{i},'_780nm_power2_stepsize1_zoom2_kalman3_Ch1.stk'];
end

for i = 1:length(on_fileCode)    
    on_fileName{i} = ['25x_',on_fileCode{i},'_power2_stepsize1_PMT700_zoom2_kalman3_Ch1.stk'];
    on_fileName1{i} = ['25x_',on_fileCode{i},'_780nm_power2_stepsize1_zoom2_kalman3_Ch1.stk'];
    if i == 22
        on_fileName1{i} = ['25x_',on_fileCode{i},'_780nm_power2_stepsize1_zoom2_Ch1.stk'];
    end
end

for i = 1:length(mix_fileCode)    
    mix_fileName{i} = ['25x_',mix_fileCode{i},'_power2_stepsize1_PMT700_zoom2_kalman3_Ch1.stk'];
    mix_fileName1{i} = ['25x_',mix_fileCode{i},'_780nm_power2_stepsize1_zoom2_kalman3_Ch1.stk'];
end

%% load Ch1 maxImages (the maximum plane)

I_ch_col = zeros(800,800,length(ch_fileCode));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_ch_col(:,:,i) = I(:,:,ind);
    figure(1);subplot(7,7,i);imagesc(I_ch_col(:,:,i),[130 2700]);axis off
end

I_on_col = zeros(800,800,length(on_fileCode));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_on_col(:,:,i) = I(:,:,ind);    
%     I_on_col(i) = mean(I(:));
    figure(2);subplot(7,7,i);imagesc(I_on_col(:,:,i),[130 2000]);axis off
end

I_mix_col = zeros(800,800,length(mix_fileCode));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_mix_col(:,:,i) = I(:,:,ind);    
%     I_mix_col(i) = mean(I(:));
    figure(3);subplot(7,7,i);imagesc(I_mix_col(:,:,i));axis off
end

%% load Ch1 max projection Images (the maximum plane)

I_ch_col = zeros(800,800,length(ch_fileCode));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    I_ch_col(:,:,i) = max(I,[],3);
    figure(1);subplot(7,7,i);imagesc(I_ch_col(:,:,i));axis off
end

I_on_col = zeros(800,800,length(on_fileCode));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end    
    I_on_col(:,:,i) = max(I,[],3);
%     I_on_col(i) = mean(I(:));
    figure(2);subplot(7,7,i);imagesc(I_on_col(:,:,i));axis off
end

I_mix_col = zeros(800,800,length(mix_fileCode));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    I_mix_col(:,:,i) = max(I,[],3);
%     I_mix_col(i) = mean(I(:));
    figure(3);subplot(7,7,i);imagesc(I_mix_col(:,:,i));axis off
end

%% quantify the Ch1 intensity

I_ch_col = zeros(length(ch_fileCode),1);
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    I_ch_col(i) = mean(I(:));
end

I_on_col = zeros(length(on_fileCode),1);
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    I_on_col(i) = mean(I(:));
end

I_mix_col = zeros(length(mix_fileCode),1);
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    I_mix_col(i) = mean(I(:));
end

%% load Ch1 maxImages (the maximum plane), calculate covariance image and intensity histogram

I_ch_col = zeros(800,800,length(ch_fileCode));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_ch_col(:,:,i) = I(:,:,ind);
%     I_ch_col(i) = mean(I(:));
end

I_on_col = zeros(800,800,length(on_fileCode));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_on_col(:,:,i) = I(:,:,ind);    
%     I_on_col(i) = mean(I(:));
end

I_mix_col = zeros(800,800,length(mix_fileCode));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_mix_col(:,:,i) = I(:,:,ind);    
%     I_mix_col(i) = mean(I(:));
end

% calculate covariance images
for i = 1:size(I_ch_col,3)
    ch_col_cov(:,:,i) = cov(double(I_ch_col(:,:,i)));
end
for i = 1:size(I_on_col,3)
    on_col_cov(:,:,i) = cov(double(I_on_col(:,:,i)));
end
for i = 1:size(I_mix_col,3)
    mix_col_cov(:,:,i) = cov(double(I_mix_col(:,:,i)));
end

save maxImage_cov ch_col_cov on_col_cov mix_col_cov
figure;
for i = 1:size(ch_col_cov,3);
    1:subplot(7,7,i);
    imagesc(ch_col_cov(:,:,i));
    axis off
end
figure;
for i = 1:size(on_col_cov,3);
    1:subplot(7,7,i);
    imagesc(on_col_cov(:,:,i));
    axis off
end

% calculate histgram in the max images
I_ch_col = reshape(I_ch_col,800*800,47);
I_on_col = reshape(I_on_col,800*800,45);
I_mix_col = reshape(I_mix_col,800*800,4);

X = 0:4095;
H_ch = zeros(length(X),size(I_ch_col,2));
for i = 1:size(I_ch_col,2)
    H_ch(:,i) = hist(I_ch_col(:,i),X);
end
H_on = zeros(length(X),size(I_on_col,2));
for i = 1:size(I_on_col,2)
    H_on(:,i) = hist(I_on_col(:,i),X); 
end
H_mix = zeros(length(X),size(I_mix_col,2));
for i = 1:size(I_mix_col,2)
    H_mix(:,i) = hist(I_mix_col(:,i),X);
end

save maxImage_hist H_ch H_on H_mix X

%% load Ch1 all image planes, find collagen volume per plane without smoothing

I_ch_col = uint16(zeros(800,800,length(ch_fileCode)));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
%     I = I - min(I(:));
    BW = zeros(size(I));    
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end
    ch_col_vol(i) = sum(BW(:))/size(BW,3);
end

I_on_col = uint16(zeros(800,800,length(on_fileCode)));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end
    end
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);      
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);        
            k1 = k1+1;
        end
    end
%     I = I - min(I(:));    
    BW = zeros(size(I));
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end    
    on_col_vol(i) = sum(BW(:))/size(BW,3);
end

I_mix_col = uint16(zeros(800,800,length(mix_fileCode)));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
%     I = I - min(I(:));    
    BW = zeros(size(I));    
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end
    mix_col_vol(i) = sum(BW(:))/size(BW,3);
end
save totalImage_col_vol ch_col_vol on_col_vol mix_col_vol

%% load Ch1 all image planes, find collagen volume per plane after smoothing

I_ch_col = uint16(zeros(800,800,length(ch_fileCode)));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    H = fspecial('gaussian',[2,2],1);
    for i = 1:size(I,3)
        I(:,:,i) = filter2(H,I(:,:,i));
    end
    BW = zeros(size(I));    
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end
    ch_col_vol(i) = sum(BW(:))/size(BW,3);
end

I_on_col = uint16(zeros(800,800,length(on_fileCode)));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end
    end
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);      
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);        
            k1 = k1+1;
        end
    end
    H = fspecial('gaussian',[2,2],1);
    for i = 1:size(I,3)
        I(:,:,i) = filter2(H,I(:,:,i));
    end
    BW = zeros(size(I));
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end    
    on_col_vol(i) = sum(BW(:))/size(BW,3);
end

I_mix_col = uint16(zeros(800,800,length(mix_fileCode)));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    H = fspecial('gaussian',[2,2],1);
    for i = 1:size(I,3)
        I(:,:,i) = filter2(H,I(:,:,i));
    end  
    BW = zeros(size(I));    
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end
    mix_col_vol(i) = sum(BW(:))/size(BW,3);
end
save totalImage_col_vol ch_col_vol on_col_vol mix_col_vol

%% load Ch1 all image planes, find collagen volume using global threshold without smoothing

I_ch_col = uint16(zeros(800,800,length(ch_fileCode)));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    planes = size(I,3);
    I = reshape(I,size(I,1),size(I,2)*size(I,3));
%     I = I - min(I(:));
    BW = zeros(size(I));
	level = graythresh(I);
	BW = im2bw(I,level);
    ch_col_vol(i) = sum(BW(:))/planes;
end

I_on_col = uint16(zeros(800,800,length(on_fileCode)));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end
    end
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);      
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);        
            k1 = k1+1;
        end
    end
    planes = size(I,3);
    I = reshape(I,size(I,1),size(I,2)*size(I,3));
%     I = I - min(I(:));
    BW = zeros(size(I));
	level = graythresh(I);
	BW = im2bw(I,level);
    on_col_vol(i) = sum(BW(:))/planes;
end

I_mix_col = uint16(zeros(800,800,length(mix_fileCode)));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    planes = size(I,3);
    I = reshape(I,size(I,1),size(I,2)*size(I,3));
%     I = I - min(I(:));
    BW = zeros(size(I));
	level = graythresh(I);
	BW = im2bw(I,level);
    mix_col_vol(i) = sum(BW(:))/planes;
end
save totalImage_col_vol ch_col_vol on_col_vol mix_col_vol

%% load Ch1 maxImages (the maximum plane), find collagen volume without smoothing

I_ch_col = uint16(zeros(800,800,length(ch_fileCode)));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_ch_col(:,:,i) = I(:,:,ind);  
end

I_on_col = uint16(zeros(800,800,length(on_fileCode)));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_on_col(:,:,i) = I(:,:,ind);  
end

I_mix_col = uint16(zeros(800,800,length(mix_fileCode)));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    [MAX,ind] = max(sum(sum(I)));
    I_mix_col(:,:,i) = I(:,:,ind);
end

for i = 1:size(I_ch_col,3)
    level = graythresh(I_ch_col(:,:,i));
	BW = im2bw(I_ch_col(:,:,i),level);
    ch_col_vol(i) = sum(BW(:));
end
for i = 1:size(I_on_col,3)
    level = graythresh(I_on_col(:,:,i));
	BW = im2bw(I_on_col(:,:,i),level);
    on_col_vol(i) = sum(BW(:));
end
for i = 1:size(I_mix_col,3)
    level = graythresh(I_mix_col(:,:,i));
	BW = im2bw(I_mix_col(:,:,i),level);
    mix_col_vol(i) = sum(BW(:));
end
save maxImage_col_vol ch_col_vol on_col_vol mix_col_vol

%% load Ch1 max projection image, find collagen volume

I_ch_col = uint16(zeros(800,800,length(ch_fileCode)));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end    
    I_ch_col(:,:,i) = max(I,[],3);
end

I_on_col = uint16(zeros(800,800,length(on_fileCode)));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    I_on_col(:,:,i) = max(I,[],3);
end

I_mix_col = uint16(zeros(800,800,length(mix_fileCode)));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = [];    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    I_mix_col(:,:,i) = max(I,[],3);
end

for i = 1:size(I_ch_col,3)
    level = graythresh(I_ch_col(:,:,i));
	BW = im2bw(I_ch_col(:,:,i),level);
    ch_col_vol(i) = sum(BW(:));
end
for i = 1:size(I_on_col,3)
    level = graythresh(I_on_col(:,:,i));
	BW = im2bw(I_on_col(:,:,i),level);
    on_col_vol(i) = sum(BW(:));
end
for i = 1:size(I_mix_col,3)
    level = graythresh(I_mix_col(:,:,i));
	BW = im2bw(I_mix_col(:,:,i),level);
    mix_col_vol(i) = sum(BW(:));
end
save maxProjImage_col_vol ch_col_vol on_col_vol mix_col_vol

%% FFT
ch_col_fft = zeros(800,800,length(ch_fileCode));
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I0 = tiffread2(filePath);         
    I1 = []; 
    k1 = 1;
    for k = 1:length(I0)
        if mean(I0(k).data(:)) > 100
            I1(:,:,k1) = uint16(I0(k).data);
            k1 = k1+1;
        end            
    end 
    I = zeros(800,800,size(I1,3));    
    for k = 1:size(I1,3)
        I(:,:,k) = log(abs(fftshift(fft2(I1(:,:,k)))));
    end 
%     I2 = mean(I,3);
%     [i0,j0]=find(I2==max(I2(:)));    
    ch_col_fft(:,:,i) = mean(I,3);
end

on_col_fft = zeros(800,800,length(on_fileCode));
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I0 = tiffread2(filePath);
    k1 = 1;
    I1 = [];
    for k = 1:length(I0)
        if mean(I0(k).data(:)) > 100
            I1(:,:,k1) = uint16(I0(k).data);
            k1 = k1+1;
        end
    end
    I = zeros(800,800,size(I1,3));
    for k = 1:size(I1,3)
        I(:,:,k) = log(abs(fftshift(fft2(I1(:,:,k)))));
    end
    on_col_fft(:,:,i) = mean(I,3);
end

mix_col_fft = zeros(800,800,length(mix_fileCode));
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I0 = tiffread2(filePath);
    k1 = 1;
    I1 = [];    
    for k = 1:length(I0)
        if mean(I0(k).data(:)) > 100
            I1(:,:,k1) = uint16(I0(k).data);
            k1 = k1+1;
        end
    end
    I = zeros(800,800,size(I1,3));
    for k = 1:size(I1,3)
        I(:,:,k) = log(abs(fftshift(fft2(I1(:,:,k)))));
    end
    mix_col_fft(:,:,i) = mean(I,3);
end
save data_fft ch_col_fft on_col_fft mix_col_fft

% plot
figure; hold on;
for i = 1:size(ch_col_fft,3)
    I = ch_col_fft(:,:,i);
    [i0,j0] = find(I==max(I(:)));
    plot(I(i0,:),'b');
end

figure; hold on;
for i = 1:size(on_col_fft,3)
    I = on_col_fft(:,:,i);
    [i0,j0] = find(I==max(I(:)));
    plot(I(i0,:),'r');
end

%% fit FFT to Lorentz
a_ch_x = zeros(size(ch_col_fft,3),1);
x0_ch_x = a_ch_x;
gamma_ch_x = a_ch_x;
a_ch_y = a_ch_x;
x0_ch_y = a_ch_x;
gamma_ch_y = a_ch_x;
for i = 1:size(ch_col_fft,3)
    I = ch_col_fft(:,:,i);
    [i0,j0] = find(I==max(I(:)));
    [a_tmp,x0_tmp,gamma_tmp]=fitlorentz(1:size(I,1),I(i0,:));
    a_ch_x(i) = a_tmp;
    x0_ch_x(i) = x0_tmp;
    gamma_ch_x(i) = gamma_tmp;
    [a_tmp,x0_tmp,gamma_tmp]=fitlorentz(1:size(I,2),I(:,j0)');
    a_ch_y(i) = a_tmp;
    x0_ch_y(i) = x0_tmp;
    gamma_ch_y(i) = gamma_tmp;        
end

a_on_x = zeros(size(on_col_fft,3),1);
x0_on_x = a_on_x;
gamma_on_x = a_on_x;
a_on_y = a_on_x;
x0_on_y = a_on_x;
gamma_on_y = a_on_x;
for i = 1:size(on_col_fft,3)
    I = on_col_fft(:,:,i);
    [i0,j0] = find(I==max(I(:)));
    [a_tmp,x0_tmp,gamma_tmp]=fitlorentz(1:size(I,1),I(i0,:));
    a_on_x(i) = a_tmp;
    x0_on_x(i) = x0_tmp;
    gamma_on_x(i) = gamma_tmp;    
    [a_tmp,x0_tmp,gamma_tmp]=fitlorentz(1:size(I,2),I(:,j0)');    
    a_on_y(i) = a_tmp;
    x0_on_y(i) = x0_tmp;
    gamma_on_y(i) = gamma_tmp;        
end

a_mix_x = zeros(size(mix_col_fft,3),1);
x0_mix_x = a_mix_x;
gamma_mix_x = a_mix_x;
a_mix_y = a_mix_x;
x0_mix_y = a_mix_x;
gamma_mix_y = a_mix_x;
for i = 1:size(mix_col_fft,3)
    I = mix_col_fft(:,:,i);
    [i0,j0] = find(I==max(I(:)));
    [a_tmp,x0_tmp,gamma_tmp]=fitlorentz(1:size(I,1),I(i0,:));
    a_mix_x(i) = a_tmp;
    x0_mix_x(i) = x0_tmp;
    gamma_mix_x(i) = gamma_tmp;    
    [a_tmp,x0_tmp,gamma_tmp]=fitlorentz(1:size(I,2),I(:,j0)');
    a_mix_y(i) = a_tmp;
    x0_mix_y(i) = x0_tmp;
    gamma_mix_y(i) = gamma_tmp;    
end
save fft_width a_ch_x x0_ch_x gamma_ch_x a_on_x x0_on_x gamma_on_x a_mix_x x0_mix_x gamma_mix_x a_ch_y x0_ch_y gamma_ch_y a_on_y x0_on_y gamma_on_y a_mix_y x0_mix_y gamma_mix_y



%% load Ch1 all image planes, find collagen area in each plane without smoothing 
%  and find perimeter using smoothed thresholded image, calculate average of area/perimeter of each plane

ch_col_shape = zeros(length(ch_fileCode),1);
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
%     I = I - min(I(:));
    BW = uint16(zeros(size(I)));
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end
    
    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    pars = struct([]);
    for k = 1:size(I,3)
        I(:,:,k) = uint16(filter2(H,I(:,:,k)));
        stats = regionprops(I(:,:,k)>0,'Area','Perimeter');
        for k1 = 1:length(stats)
            Area(k,k1) = stats(k1).Area;
            Perimeter(k,k1) = stats(k1).Perimeter;
        end        
    end

    Area_ch(i) = sum(Area(:));
    Perimeter_ch(i) = sum(Perimeter(:));
    
%     ch_col_shape(i) = Area_ch(i)./Perimeter_ch(i);    
    ch_col_shape(i) = squeeze(sum(sum(sum(BW))))./Perimeter_ch(i);
end

on_col_shape = zeros(length(on_fileCode),1);
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end
    end
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);      
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);        
            k1 = k1+1;
        end
    end
%     I = I - min(I(:));    
    BW = zeros(size(I));
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end    

    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    Area = [];
    Perimeter = [];
    for k = 1:size(I,3)
        I(:,:,k) = uint16(filter2(H,I(:,:,k)));
        stats = regionprops(I(:,:,k)>0,'Area','Perimeter');
        for k1 = 1:length(stats)
            Area(k,k1) = stats(k1).Area;
            Perimeter(k,k1) = stats(k1).Perimeter;
        end        
    end

    Area_on(i) = sum(Area(:));
    Perimeter_on(i) = sum(Perimeter(:));
    
%     on_col_shape(i) = Area_on(i)./Perimeter_on(i);    
    on_col_shape(i) = squeeze(sum(sum(sum(BW))))./Perimeter_on(i);
end

mix_col_shape = zeros(length(mix_fileCode),1);
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
%     I = I - min(I(:));    
    BW = zeros(size(I));    
    for k = 1:size(I,3)
        level = graythresh(I(:,:,k));
        BW(:,:,k) = im2bw(I(:,:,k),level);
    end

    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    Area = [];
    Perimeter = [];
    for k = 1:size(I,3)
        I(:,:,k) = uint16(filter2(H,I(:,:,k)));
        stats = regionprops(I(:,:,k)>0,'Area','Perimeter');
        for k1 = 1:length(stats)
            Area(k,k1) = stats(k1).Area;
            Perimeter(k,k1) = stats(k1).Perimeter;
        end        
    end

    Area_mix(i) = sum(Area(:));
    Perimeter_mix(i) = sum(Perimeter(:));
    
%     mix_col_shape(i) = Area_mix(i)./Perimeter_mix(i);    
    mix_col_shape(i) = squeeze(sum(sum(sum(BW))))./Perimeter_mix(i);
end
save totalImage_col_shape ch_col_shape on_col_shape mix_col_shape

%% load Ch1 max image plane, find collagen area without smoothing 
%  and find perimeter using smoothed thresholded image, calculate area/perimeter

ch_col_shape = zeros(length(ch_fileCode),1);
for i = 1:length(ch_fileCode)
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end
    end
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    
    [MAX,ind] = max(sum(sum(I)));
    I = I(:,:,ind);    

    BW = uint16(zeros(size(I)));
    level = graythresh(I);
    BW = im2bw(I,level);
    
    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    
    I = uint16(filter2(H,I));
    stats = regionprops(I>0,'Area','Perimeter');
    
    for k = 1:length(stats)
        Area(k) = stats(k).Area;
        Perimeter(k) = stats(k).Perimeter;
    end
    Area_ch(i) = sum(Area);
    Perimeter_ch(i) = sum(Perimeter);
    
%     boundaries = bwboundaries(I);
%     B = [];
%     for k = 1:length(boundaries)
%         B(k) = length(boundaries{k});
%     end
%     Perimeter_ch(i) = sum(B);
    ch_col_shape(i) = mean(Area_ch(i)./Perimeter_ch(i));
%     ch_col_shape(i) = squeeze(sum(sum(BW)))./Perimeter_ch(i);
%     figure(3);subplot(7,7,i);imagesc(I);axis off
end

on_col_shape = zeros(length(on_fileCode),1);
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end
    end
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);      
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);        
            k1 = k1+1;
        end
    end
    
    [MAX,ind] = max(sum(sum(I)));
    I = I(:,:,ind);    

    BW = uint16(zeros(size(I)));
    level = graythresh(I);
    BW = im2bw(I,level);
    
    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    pars = struct([]);
    I = uint16(filter2(H,I));
    stats = regionprops(I>0,'Area','Perimeter');
    
    for k = 1:length(stats)
        Area(k) = stats(k).Area;
        Perimeter(k) = stats(k).Perimeter;
    end
    Area_on(i) = sum(Area);
    Perimeter_on(i) = sum(Perimeter);
    
%     boundaries = bwboundaries(I);
%     B = [];
%     for k = 1:length(boundaries)
%         B(k) = length(boundaries{k});
%     end
%     Perimeter_on(i) = sum(B);    
	on_col_shape(i) = Area_on(i)./Perimeter_on(i);
% 	on_col_shape(i) = squeeze(sum(sum(BW)))./Perimeter_on(i);
%     figure(3); subplot(7,7,i); imagesc(BW); axis off
%     figure(4); subplot(7,7,i); imagesc(I>0); axis off    
%     figure(6);subplot(7,7,i);imagesc(I);axis off
end

mix_col_shape = zeros(length(mix_fileCode),1);
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end

    [MAX,ind] = max(sum(sum(I)));
    I = I(:,:,ind);    

    BW = uint16(zeros(size(I)));
    level = graythresh(I);
    BW = im2bw(I,level);
    
    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    pars = struct([]);
    I = uint16(filter2(H,I));
    stats = regionprops(I>0,'Area','Perimeter');
    
    for k = 1:length(stats)
        Area(k) = stats(k).Area;
        Perimeter(k) = stats(k).Perimeter;
    end
    Area_mix(i) = sum(Area);
    Perimeter_mix(i) = sum(Perimeter);
    
%     boundaries = bwboundaries(I);
%     B = [];
%     for k = 1:length(boundaries)
%         B(k) = length(boundaries{k});
%     end
%     Perimeter_mix(i) = sum(B);    
    mix_col_shape(i) = Area_mix(i)./Perimeter_mix(i);
%     mix_col_shape(i) = squeeze(sum(sum(BW)))./Perimeter_mix(i);
end
save totalImage_col_shape ch_col_shape on_col_shape mix_col_shape


%% load Ch1 max Projection image, find collagen area without smoothing 
%  and find perimeter using smoothed thresholded image, calculate area/perimeter

ch_col_shape = zeros(length(ch_fileCode),1);
for i = 1:length(ch_fileCode)        
    i
    for j = 1:length(ch_fileNumber)
        if sum(ch_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end    
    filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' ch_folder{folderNo} '\Processed images\tumor\' ch_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);    
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end
    
    I = max(I,[],3);
    
    BW = uint16(zeros(size(I)));
    level = graythresh(I);
    BW = im2bw(I,level);
    
    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    pars = struct([]);
    I = uint16(filter2(H,I));
    stats = regionprops(I>0,'Area','Perimeter');

    for k = 1:length(stats)
        Area(k) = stats(k).Area;
        Perimeter(k) = stats(k).Perimeter;
    end
    Area_ch(i) = sum(Area);
    Perimeter_ch(i) = sum(Perimeter);    
%     ch_col_shape(i) = mean(Area_ch(i)./Perimeter_ch(i));    
    ch_col_shape(i) = mean(squeeze(sum(sum(BW)))./Perimeter_ch(i));    
%     figure(3);subplot(7,7,i);imagesc(I);axis off
end

on_col_shape = zeros(length(on_fileCode),1);
for i = 1:length(on_fileCode)
    i
    for j = 1:length(on_fileNumber)
        if sum(on_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end
    end
    filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' on_folder{folderNo} '\Processed images\tumor\' on_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);      
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);        
            k1 = k1+1;
        end
    end
    
    I = max(I,[],3);

    BW = uint16(zeros(size(I)));
    level = graythresh(I);
    BW = im2bw(I,level);
    
    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    pars = struct([]);
    I = uint16(filter2(H,I));
    stats = regionprops(I>0,'Area','Perimeter');

    for k = 1:length(stats)
        Area(k) = stats(k).Area;
        Perimeter(k) = stats(k).Perimeter;
    end
    Area_on(i) = sum(Area);
    Perimeter_on(i) = sum(Perimeter);    
% 	on_col_shape(i) = mean(Area_on(i)./Perimeter_on(i));
	on_col_shape(i) = mean(squeeze(sum(sum(BW)))./Perimeter_on(i));
%     figure(3); subplot(7,7,i); imagesc(BW); axis off
%     figure(4); subplot(7,7,i); imagesc(I>0); axis off    
%     figure(6);subplot(7,7,i);imagesc(I);axis off
end

mix_col_shape = zeros(length(mix_fileCode),1);
for i = 1:length(mix_fileCode)
    i
    for j = 1:length(mix_fileNumber)
        if sum(mix_fileNumber(1:j)) >= i
            folderNo = j;
            break
        end        
    end        
    filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName{i}];
    if ~exist(filePath)
        filePath = ['Onco-chromo\' mix_folder{folderNo} '\Processed images\tumor\' mix_fileName1{i}];
    end
    I1 = tiffread2(filePath);
    k1 = 1;
    I = uint16([]);
    for k = 1:length(I1)
        if mean(I1(k).data(:)) > 100
            I(:,:,k1) = uint16(I1(k).data);
            k1 = k1+1;
        end
    end

    I = max(I,[],3);

    BW = uint16(zeros(size(I)));
    level = graythresh(I);
    BW = im2bw(I,level);
    
    I = I.*uint16(BW);
    H = fspecial('gaussian',[10,10],3);
    pars = struct([]);
    I = uint16(filter2(H,I));
    stats = regionprops(I>0,'Area','Perimeter');

    for k = 1:length(stats)
        Area(k) = stats(k).Area;
        Perimeter(k) = stats(k).Perimeter;
    end
    Area_mix(i) = sum(Area);
    Perimeter_mix(i) = sum(Perimeter);    
%     mix_col_shape(i) = mean(Area_mix(i)./Perimeter_mix(i));
    mix_col_shape(i) = mean(squeeze(sum(sum(BW)))./Perimeter_ch(i));
end
save totalImage_col_shape ch_col_shape on_col_shape mix_col_shape


% 
% boundaries = bwboundaries(Ibw_smooth(:,:,1)>0);
% hold on;
% for k=1:length(boundaries)
%    b = boundaries{k};
%    plot(b(:,2),b(:,1),'b','LineWidth',1);
% end

% --------------------------------- lorentz.m ---------------------------------------%
function y = lorentz(x,x0,gamma)
y = 1/pi*0.5*gamma./((x-x0).^2+(0.5*gamma)^2);
plot(x,y);
drawnow

% --------------------------------- fitlorentz.m ---------------------------------------%
function [a,x0,gamma] = fitlorentz(x,y)
err = @(pars,x,y) sqrt(sum((y - pars(1)*lorentz(x,pars(2),pars(3))-pars(4)).^2));
[MAX,ind] = max(y);
a = sum(y(:));
x0 = x(ind);
gamma = (x(end)-x(1))/10;
y0 = min(y);
pars = fminsearch(@(pars)err(pars,x,y),[a,x0,gamma,y0]);
a = pars(1)
x0 = pars(2);
gamma = pars(3);
y0 = pars(4);
yfit = lorentz(x,x0,gamma);
plot(x,y,'r.',x,a*yfit+y0)
drawnow

% --------------------------------- gaussian.m ---------------------------------------%
function y = gaussian(x,a,mu,sigma)
y = a*exp(-(x-mu).^2/2/sigma^2);

% --------------------------------- fitgaussian.m ---------------------------------------%
function [mu,sigma,yfit] = fitgaussian(x,y)
err = @(pars,x,y) sqrt(sum((y - gaussian(x,pars(1),pars(2),pars(3))).^2));
[MAX,ind] = max(y);
pars = fminsearch(@(pars)err(pars,x,y),[1,x(ind),x(2)-x(1)]);
a = pars(1)
mu = pars(2);
sigma = pars(3);
yfit = gaussian(x,a,mu,sigma);
bar(x,y)
hold on
plot(x,yfit)